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frequencies  of  radially  symmetric  modes  of  vibration  are 
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spherical  cavity.  To  minimize  experimental  perturbations  a 

relative  measurement  against  argon  gas  used  as  a  standard 

reference  was  employed.  The  sonic  speed  in  turn  is  utilized  as 

a  sensitive  probe  for  detecting  phase  boundaries  and  critical 

points  of  substances.  Location  of  the  boundary  is  revealed  by 

a  discontinuity   in   the  variation   of   sonic   speed  with 

temperature  for  a  sample  confined  isochorically  in  the 

resonator.   The  corresponding   temperature   at  which   this 

discontinuity   occurs   was   taken   as   a   phase   boundary 

temperature.  This  indicator  is  very  pronounced,  especially 
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different  loading  densities  around  the  critical  value. 
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CHAPTER  1 
INTRODUCTION 


Natural  gas  and  gasoline  are  primarily  mixtures  of  the 
lighter  hydrocarbons  with  varying  amounts  of  nonhydrocarbons 
such  as  water,  carbon  dioxide,  and  hydrogen  disulfide.  Heavier 
fossil  fuel  mixtures  such  as  crude  oil  consist  of  a  myriad  of 
higher  boiling  hydrocarbons  and  various  compounds  containing 
sulfur,  nitrogen,  and  oxygen.  In  designing  production, 
processing,  transport,  and  handling  systems  for  these 
materials  reliable  knowledge  of  their  physical  properties  is 
crucial.  For  example,  in  the  important  activities  of  enhanced 
petroleum  recovery  by  carbon  dioxide  injection  into  a 
resevoir1  and  transport  of  natural  gas  mixtures  containing 
carbon  dioxide,  precise  knowledge  of  the  equilibrium  phase 
behaviors  of  the  systems  is  important.  Additionally,  the  need 
exists  for  knowledge  of  critical  point  parameters,  hydrate 
formation  conditions,  density,  enthalpy,  dew  points,  bubble 
points,  etc.  Information  of  these  kinds  about  select  pure 


1Under  the  right  condition  of  pressure,  temperature,  and 
oil  composition,  carbon  dioxide  will  mix  with  the  crude  in  the 
reservoir  to  form  a  single  phase  liquid  which  is  much  lighter 
than  the  original  oil  and  consequently  easier  to  bring  to  the 
surface. 
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components  and  mixtures  is  further  needed  to  support  the 

development  of  semiempirical  equations  of  state  and  other 

correlations  useful  for  predicting  the  behavior  of  a  wide 

range  of  systems.  For  this  reason,  much  research  has  been  done 

in  both  the  measurement  (experimental  method)  and  prediction 

(theoretical  method)  of  hydrocarbon  fluid  properties.  The 

great  variety  of  chemical  systems  of  practical  interest, 

combined    with   the   experimental   difficulties   and   time 

investment  involved  with  direct  measurements,  is  a  powerful 

driving  force  behind  the  development  of  useful  correlations. 

Numerous  empirical  equations  of  state  have  been  developed 

over  the  years.  The  pioneering  work  of  van  der  Waals  focused 

attention  on  the  possible  causes  for  nonideal  gas  behavior. 

Although  his  efforts  remain  historically  instructive,  his 

equation  of   state(1)   is   little  used  today  having  been 

superseded  by  a  succession  of  others  ranging  in  complexity 

from  those  of  Benedict,  Webb  and  Rubin, (2)  Redlich  and  Kwong,(3) 

and  Peng  and  Robinson'4'  to  much  more  sophisticated  ones 

containing  forty  or  more  coefficients  such  as  the  AGA8 

equation, <5)  for  example.  Also  a  number  of  semi-theoretical  and 

theoretical  equations  based  on  molecular  models  have  been 

devised  (see  for  examples  references  6-11.)  They  were  invented 

to  rectify  for  non-ideality  both  pure  and  multicomponent  gases 

and,  in  general,  they  work  remarkably  well  for  predictions  of 

gas  properties  in  the  noncritical  region.  They  generally  fail, 

however,  to  describe  the  unique  behavior  of  fluids  in  the 
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critical  region,  but  progress  toward  overcoming  this 
shortcoming  is  being  made  on  some  fronts. (12_15>  Thus  far  all 
models  developed  to  work  in  the  critical  region  require 
knowledge  of  the  critical  parameters;  therefore  the  prediction 
abilities  of  these  models  rely  partly  on  the  accuracy  of  these 
parameters  which  in  turn  must  be  obtained  experimentally.  This 
is  typically  done  only  for  pure  components,  and  the  critical 
behavior  of  mixtures  is  deduced  using  various  combining 
schemes.  The  chance  for  introducing  uncertainty  increases  then 
as  the  complexity  of  the  mixture  and  is  generally  substantial 
for  most  mixtures  of  importance  to  the  gas  industry.  One 
popular  combining  scheme  uses  Kay's  mixing  rules*16':  P  =EY.P  . 
and  T  =SY.T  .  where  P   and  T   are  the  so-called  pseudocritical 

pc       1   CI  pc  pc  m 

pressure  and  temperature,  respectively,  for  the  mixture,  and 
Y±,  Pci,  Tci  is  the  mole  fraction,  pressure,  and  temperature, 
respectively,  of  component  i  in  the  mixture.  Alternatively, 
pseudocritical  parameters  may  be  obtained  from  a  correlation 
based  on  a  collection  of  natural  gas  data'17'18'.  The  values  for 
a  specific  mixture  in  general  depend  on  the  combining  scheme 
used  to  predict  them. 

The  importance  of  accurate  values  of  the  critical  point 
parameters  of  both  pure  and  mixed  fluids  is  established.  Good 
values  for  many  pure  components  are  available,  but  very  few 
mixtures  have  received  careful  experimental  study.  The  higher 
the  number  of  system  components  the  fewer  the  number  of 
tabulated  values  are  generally  available.  This  scarcity  is 


4 
partly  due  to  the  difficulty  in  location  of  the  critical  point 
in  mixtures  by  conventional  methods.  The  most  widely  used 
method  is  based  on  direct  observation  of  the  appearance  and 
disappearance  of  a  meniscus,  indicating  a  boundary  line 
between  gas  and  liquid  phases  when  a  system  enters  the 
critical  condition. (19)  This  method  relies  on  the  visual  acuity 
of  the  operator,  and  clearly  has  some  subjectivity. 

Acoustic  determination  of  thermodynamic  properties  such 
as  heat  capacity,  virial  coefficients,  vapor  pressure,  etc.  by 
spherical  resonator  techniques  has  been  very  successful  in  the 
past,  especially  in  this  research  group. (20_26)  This  success 
stems  mainly  from  the  unusual  sensitivity  of  the  cavity 
resonator.  Resonance  frequencies  can  be  readily  measured  to 
precisions  of  one  part  in  106  or  better  compared  to  one  part 
in  104  or  105  for  pressure,  temperature,  volume  density  etc. 
The  acoustic  resonator  also  benefits  thermodynamic 
measurements  by  operating  at  low  frequencies.  At  high 
frequencies  the  speed  of  sound  is  affected  by 
irreversibilities  resulting  principally  from  the  delay  of 
energy  flow  into  and  out  of  internal  energy  modes, 
particularly  vibrational  modes. (27) 

This  work  is  a  first  attempt  to  bring  the  superiority  of 
the  spherical  acoustic  resonance  technique  to  bear  on 
development  of  a  new  approach  to  critical  point  determination. 
Experiments  were  performed  according  to  logical  scientific 
approach.  To  start  with,  verification  of  this  new  technique 


5 
was  made  using  a  gas  system  for  which  the  critical  point 
parameters  are  well-established.  That  gas  was  pure  carbon 
dioxide,  which  has  the  richest  history  of  prior  critical 
behavior  studies.  Choosing  this  gas  is  beneficial  not  only 
because  of  its  well-known  values  of  critical  point  parameters 
but  also  for  its  recognized  importance  in  the  fuel  industry  as 
a  near-critical-state  solvent  in  enhanced  oil  recovery  and  as 
a  supercritical  extractant.  Next,  the  complexity  of  the  gas 
system  was  increased  by  selecting  a  binary  component  gas 
mixture  (C02~C2H6)  for  study.  The  lessons  learned  in  working 
with  these  first  two  systems  were  applied  to  a  follow-up  study 
of  the  other  component  of  the  mixture,  ethane,  to  further 
confirm  the  reliability  of  the  technique.  Chapter  two  will 
briefly  review  some  relevant  equations  of  acoustics  and 
explain  how  to  identify  individual  resonance  frequencies, 
which  are  the  key  parameters  leading  to  values  of  the  speed  of 
sound.  Chapter  three  will  review  some  aspects  of  the  phase 
behavior  of  both  pure  and  mixed  gases.  The  critical  point 
determination  will  subsequently  be  discussed,  starting  from  a 
background  on  conventional  approaches  to  the  newly  proposed 
approach,  which  utilizes  speed  of  sound  as  a  critical  point 
indicator.  Chapter  four  explains  comprehensively  the 
experimental  procedures  used  in  this  work.  Chapter  five 
contains  results  including  discussions  and  observations.  The 
final  chapter  is  a  conclusion. 


CHAPTER  2 
THEORETICAL  BACKGROUND 


Acoustics 


This  derivation  follows  that  given  by  references  28,  29, 
30,  and  31.  The  propagation  of  pressure  waves  in  a  lossless 
fluid  at  rest  satisfy  the  general  wave  equation'28' 

V2^  -  ±§1*  (2.1) 

C2dt2 

where  V2  =  Laplacian  operator,    <p  is  the  velocity  potential  and 
C   is   the  wave  velocity    (speed   of    sound). 
Assuming  that   <f>  =  <f>oeiot   equation    (2.1)    becomes 

Vcp  +  —<f>  -  0  (2.2) 

C2 

or 

V2^  +  k20  -  0  (2-3) 

where  k  =  u/C  =  wave  number.  If  we  consider  only  the  time- 
independent  wave,  equation  (2.3)  becomes 

V2^  +  k2<f>0   -  0  (2.4) 


If  the  variables  are  separable  0o  may  be  written  as  <p      = 
f(r)g(6)h(q>),  and  equation  (2.4)  becomes :<29) 


:2££  +  2rfJ!  +  (k2r2-M2)f  -0 
dr2     dr 

(  \ 


*  4(sin6^1  + 
¥IHB  d5     dB1 


sin26j 


g  -0 


d2h 
dq>2 


+  q2h  -  0. 


(2.5) 


where  M2  and  q2  are  separation  constants.  The  first  of  these 
three  equations  gives  the  radial  part  and  the  other  two  the 
angular  part  of  the  solution. 
Angular  Part 

The  last  equation  in  equations  (2.5)  is  a  linear 
homogeneous  second-order  differential  equation.  Its  solutions 
are  h  =  sine(q(p)  and  cosine  (q<p).  In  order  that  h  be  single- 
valued  [h(r,0,<p+2)T)  =  h(r,6,q>)]  values  of  q  must  be  integers, 
i.e.  q  =  0,  ±1,  ±2,....  Denoting  this  integer  by  m  the  second 
of  equations  (2.5)  becomes 


iTnT5dl5 


sin6 


dg 


M2  - 


m' 


sin26 


g  -  0 


(2.6) 


This  equation  has  finite  solutions  for  0  =  0  &  n   when 

M2  =  A(£  +  l)    ;  0.   =0,1,2,  ... 
Then 


1       ^  /   ±   0dg\ 
sIHedm  H0 


i(je+i)  - 


mz 


sin2e, 


g  -  0. 


(2.7) 


The  general  solution  of  this  equation  is 
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0-A,P"(cos0)  +BiQ"(cos6).  (2.8) 

where  P^fcosS)  =    The  associated  Legendre  polynomial  of  the 

first  kind  of  degree  £   and  order  m 
Q^fcose)  =    The  associated  Legendre  polynomial  of  the 
second  kind  of  degree  £   and  order  m. 
Since  Q^cosO)  becomes  infinite  at  0  =  0  &  n,    it  does  not 
apply  to  this  physical  situation'28' .  Consequently, 

0  -  A.P^cosO) 

dmP  fcosft)  (2.9) 

where      P™(cos8)   -   (sin6)m/2        *(  ' 

d(cos8)m 

if  m  =  0,  P^°(cos0)  is  the  Legendre  polynomial. 

Finally,  the  general  solution  of  the  angular  part  can  be 

written  as  follows: 

Y"(0,<p)  -P™(cos8)=°>(p).  (2.10) 

These  are  the  spherical  harmonics  of  the  first  kind. 

Radial  Part 

Equation  2.5a  is: 


-2U  A  ^  o^UJ-  ^  >  1,2^.2 


—  +  24^ 
dr2     ^ 


+  2rr±  +  (k2r2  -  M2)f  -  0.         (2.11) 


Using  the  above  result  for  the  angular  part: 

M2  =  £(£+1)    ;  £   =  0,1,2,...  equation  (2.11)  becomes 


:  +  2r^  +  (k2r2  -  £  (je  +  1)  )  f  -  0.       (2.12) 


dr 


2 


dT 


Dividing   both   sides   by   k2r2   gives 


1  d2f    2  df 


L   _   l(l  +  l)lf 


(2.13) 


k2dr2   k2r^   I,     k2r2 
Changing  the  form  of  the  function  from  f(r)  to  F(Z)  where  Z  = 
kr  and  rearranging  yields: 


*Sja  ♦  $-g«l  ♦.(!  -  ii^),,.,  .  ..     (2.14, 

Particular  solutions  are  the  spherical  Bessel  function  of 

three  kinds: 

First  kind, 

F(z)  -  3,(2)  - 


N 


JL  J  i  (Z). 

ZZ     (*+-) 


(2.15) 


Second  kind, 


F(z)  -  y,(z) 


4^-  Y   ,  (Z) 


(2.16) 


Third  kind, 

F(z)  -  h^'fz)  -  j,(z)  +  iYi(z) 


N 


4   *»,  (Z) 

ZZ     (*  +  -) 
2 


.(2) 


F(z)  -  h^'(z)  -  j,(z)  -  iYi(z)  - 


(2.17) 


JL   B<'\  (Z) 

\  ZZ     (i  +  -) 


where  J(Z),  Y(Z),  and  H(Z)  is  Bessel  function  of  the  first, 
second,  and  third  kind,  respectively. 
With,  for  example, 


Do(z) 
Ji(z) 

Yo  " 
Yi  ~ 


sin( z ) 


[sin(z)  -  z*cos(z)] 

z1 
_  cos ( z ) 

z 
[cos(z)  +  z- sin(z) ] 


(2.18) 
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For  higher  orders  of  1  the  solutions  may  be  obtained  from  the 
following  recurrence  relation: 

*|f  ♦!!(*)  "  kW    +  l)f|(2)  "  *(,-!>•  (2'19) 

The  general  solution  may  be  written  as  a  linear  combination  of 
the  first  and  second  kind  solutions  as  follows: 

F(kr)  -Aj,(kr)  +  By,(kr).  (2.20) 

F,  however,  must  be  finite  at  r=0 .  Spherical  Bessel  functions 

of  the  second  kind  do  not  satisfy  this  constraint,  so  the 

equation  (2.20)  becomes 

F(kr)  =  Aj^(kr).  (2.21) 

A  second  condition  imposes  that  the  normal  velocity  must 

vanish  at  the  rigid  wall  of  the  sphere,  i.e. 

agjkr>)   -0.  (2.22) 

This  is  a  transcendental  equation  whose  roots  give  the  normal 
frequencies  of  vibration,  £.  For  example  if 

j0(kr)  -  B^El 


,  dj.(kr) 
I  JoV ' 


sin(ka)    cos(ka)  _  Q 

ka2        a  (2.23) 

sin(ka)  _  cos(ka) 

ka2        * 
Hence,    tan(ka)  -  ka  -  ^  n. 

For  each  value  of  £  there  exists  an  infinite  number  of  $ 
values  designated  as  ^  n.  The  first  few  points  corresponding 
to  roots  of  equation  (2.23)  for  the  first  six  orders  are  shown 
in  figure  2.1. 
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Overall  Solutions 

Combining  the  results  of  the  radial  and  angular  parts 
yields  the  solution  of  the  time-independent  monochromatic  wave 
equation  in  the  spherical  polar  coordinate  system  as 

oo    +g  oo 

<Mr,6,(p)  ■JEEj.l^.rl^f).       (2.24) 

n-0  m—  t  1-0 


It  is  clear  that  in  order  to  determine  fully  any  particular 
mode  of  vibration,  we  must  specify  the  values  of  three 
characteristic  numbers:  n,  2,    and  m. 

The  roots  of  the  radial  part  £  n  correspond  to  the 
frequency  through  the  relation 

%t  ,  -  ka  -  ™   -  2nf*'"a 

r      c  C      C  (2.25) 

f        *l,n" 

where  C  =  speed  of  sound 

a  =  the  radius  of  the  spherical  cavity 
fj   =  resonance  frequencies  corresponding  to  £ 
£   =   root  of  the  spherical  Bessel  function  of  the  first 
kind 
=  eigenvalues. 
Lord  Rayleigh(30)  was  the  first  to  study  this  problem.  In 
1872  he  solved  for  the  lowest  23  eigenvalues  (sometimes  called 
natural  frequencies).  In  1952  H.  G.  Ferris<31>  revised  and 
extended  the  list  to  the  lowest  84  eigenvalues  as  shown  in 
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Fxgure  2.1.  The  first  six  orders  of  the  spherical  Bessel 
function.  Points  where  slope  equals  zero  yield  eigen  values, 
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table  2.1.  Even  more  complete  tabulations  can  be  found  in  the 
literature*32'. 


Searching  for  Resonance  Frequencies 

£  C 
From  the  relationship  f „  -  *'"    it  is  obvious  that  in 

order  to  find  the  resonance  frequency  corresponding  to  the 

normal  mode  of  vibration  of  an  eigenvalue,  ?  n,  the  speed  of 

sound  and  the  radius  of  the  sphere  must  be  known.  The  latter 

may  be  obtained  precisely  from  geometric  considerations  on  the 

resonator  cavity.  The  former  may  be  approximately  calculated 

from  an  assumed  equation  of  state.  In  this  study  we  used  the 

American  Gas  Association  equation  called  AGA8(5)  which  requires 

as  input  only  temperature,  pressure,  and  composition.  Using 

the  predicted  sonic  speed  together  with  the  measured  cavity 

radius   and  the  tabulated  eigenvalue,   the  corresponding 

resonance  frequency  ±1  n  may  be  approximately  calculated  and 

then  experimentally  searched  by  scanning  frequency  around  that 

value.  When  excited  at  a  resonance  frequency  the  amplitude  of 

the  detected  acoustic  signal  sharply  increases.  The  series  of 

resonance  peaks  arranged  by  frequency  defines  a  zeroth  order 

acoustic  spectrum  of  resonances  of  gas  enclosed  in  the 

spherical  cavity. 

Once  a  resonance  frequency  has  been  found,  it  can  serve 

as  an  internal  reference  point  to  search  for  another,  as 

follows: 
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C 


fi,n  _  C  (2.26) 

-  g  • 


If  g  is  assumed  to  be  a  constant,  this  enables  us  to  calculate 
a  resonance  frequency  corresponding  to  each  root  E        from  the 


following  relationship: 
(f,,.)r.f   f, 


"ref 


(S*,n)ref     Vref 

-  g 

(f,,n) 

at  any  other 

Consequently, 

f*,n  -  (5if„)g. 

(2.27) 


(2.28) 

Note  that  this  method  also  serves  as  a  self  check  routine  for 
the  identity  of  the  first  reference  frequency.  If  the  observed 
frequency  is  not  correctly  identified,  it  is  unlikely  that  we 
will  find  another  resonance  frequency  in  the  expected  region 
by  using  the  above  equation. 

So  far  we  have  discussed  a  zeroth  order  perturbation 
relationship: 

f,,n  "  5,,nJL  •  (2.29) 

In  reality  there  are  higher  perturbations  to  this  relation. 
Moldover  has  suggested  that  these  effects  should  be  added  to 
the  above  equation  in  the  following  manner'"'34' : 


15 
f    n  -    (^;n)(C)    +  AfB  +  Afel  +  Aft  +  Af    om...  (2.30) 

2,n  -y  ~—  a        el        t        geom 

The  first  correction  term  is  due  to  the  thermal  boundary  layer 
effect  at  the  wall  of  the  resonator.  The  second  correction 
term  is  due  to  the  finite  elastic  compliance  of  the  sphere 
which  is  not  infinitely  rigid.  The  third  term  corrects  for 
effects  of  the  gas  entrance  and  gas  exit  tubes  which 
contribute  departures  from  sphericity.  The  last  term  corrects 
for  imperfection  of  the  cavity  geometry.  For  absolute 
measurements  the  exact  treatment  could  be  computationally  and 
experimentally  demanding. 

These  problems,  however,  can  be  greatly  reduced  by 
utilizing  a  relative  measurement  technigue  as  follows: 
equation  (2.30)  may  be  written  as 

f   .  (y<c)  (2.31) 

where  \l  is  an  effective  eigenvalue  which  depends  on  the 
unperturbed  eigenvalue  and  all  perturbations  for  the 
particular  mode  (l,n).  If  it  is  assumed  to  be  independent  of 
gas  properties,  then 

.  (2ira)(fItn) 

i,n  c  V        / 

Consequently,  if  we  choose  some  system  whose  speed  of  sound  is 
known  or  can  be  calculated  precisely  from  a  reliable  equation 
of  state,  we  will  get  a  relationship: 
or 
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H%1  -  2"a(%-l 


(2.33) 


C,  "  cA 


if0 

r  T~ 


(2.34) 


where  I  refers  to  the  system  under  investigation.  The  speed  of 
sound  of  argon  is  calculated  from  the  virial  equation 
truncated  after  the  third  term  as  follows'35': 


,2    yRT 


Ar 


M. 


V     v2 


(2.35) 


where  y  =  heat  capacity  ratio  (C  /C  ) 

MAr  =  Molecular  mass  of  argon 

R  =  Gas  constant 

T  =  Absolute  Temperature 

v  =  molar  volume 

A2(T)  =  Second  acoustic  virial  coefficient 

A2(T)  =  Third  acoustic  virial  coefficient. 
The  effective  eigenvalue  determined  this  way  is  not  exact,  but 
the  effect  on  sonic  speed  measurements  is  small  and  of  no 
consequence   for   the   accurate   determination   of   phase 
boundaries. 

The  practical  advantages  of  spherical  resonators  is 
realized  for  the  purely  radial  modes  of  vibration,  which 
involve  no  tangential  motion  of  the  gas  with  respect  to  the 
resonator  wall,  i.e.  no  viscous  drag  effect.  In  addition,  the 
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radial  mode  resonance  peaks  typically  have  narrow  half -width, 
(high  Q  value2).  Furthermore,  resonance  frequencies  of  the 
radial  modes  are  sensitive  to  imperfections  of  non-sphericity 
only  to  second  order(36).  Consequently,  from  all  points  of  view 
the  most  accurate  measurements  of  the  speed  of  sound  should  be 
obtained  by  utilizing  the  radial  modes  of  vibration.  The  first 
few  radial  modes  are  mostly  used  because  of  their  associated 
low  frequencies.  A  more  complete  account  of  these  effects  can 
be  found  in  references  33  and  34. 


2 


Q  is  usually  defined  as: 

q  _  Center  Frequency  of  Resonance 
Frequency  Width  -  3dB  Points 
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CHAPTER  3 
APPLICATION  OF  SPEED  OF  SOUND  TO  LOCATE  CRITICAL  POINTS 


Phase  Behavior  of  Fluids  and  Fluid  Mixtures 

The  critical  point  of  a  fluid,  whether  it  be  a  pure 
substance  or  a  mixture,  is  a  property  of  considerable 
practical  as  well  as  theoretical  importance.  This  is  because 
the  critical  point  identifies  the  temperature,  pressure  and 
density  at  which  the  liquid  and  gaseous  phases  have  identical 
properties  and  is,  therefore,  a  key  point  in  the  construction 
of  the  phase  diagram.  Also,  a  knowledge  of  critical 
temperature  and  critical  pressure  makes  possible,  through  the 
principle  of  corresponding  states  and  an  appropriate  equation 
of  state,  the  prediction  of  the  thermodynamic  properties  of 
the  compound  when  these  properties  have  not  yet  been 
experimentally  determined.  The  need  for  information  about 
phase  behavior  has  substantially  increased  in  recent  years. 
For  example,  the  petroleum  industry  has  given  much  attention 
to  the  phase  behavior  of  fluid  mixtures  at  high  pressures  to 
support  useful  strategies  to  increase  the  yield  of  petroleum 
from  oil  wells.  To  have  a  clear  picture  of  the  phase  behavior 
at  high  pressure  and  temperature  of  a  fluid  mixture,  the 
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critical  point,  the  maximum  pressure  and  the  maximum 
temperature  for  the  heterogeneous  region  must  be  known. 

Evidence  for  the  existence  of  a  critical  point  was  first 
presented  in  1823  by  de  La  Tour<37),  who  observed  that  a  liquid 
when  heated  in  a  hermetically  sealed  glass  tube  is  reduced  to 
vapor  in  a  volume  from  two  to  four  times  the  original  volume 
of  the  sample.  However,  it  was  not  until  the  quantitative 
measurements  of  Andrews  on  carbon  dioxide'38'39'  in  1869  that  the 
nature  of  the  transition  was  understood.  He  was  the  first  to 
coin  the  term  "critical  point"  for  the  phenomenon  associated 
with  this  liquid-vapor  transition.  For  fluid  mixtures  the 
first  reliable  experimental  investigation  of  the  critical 
state  began  with  the  work  of  Kuenen'40'41'  in  1897.  Interest  in 
the  critical  region  in  the  period  1876-1914  was  peaked  by  the 
experimental  and  theoretical  studies  of  pressure-volume- 
temperature  (p-v-T)  relationship  for  both  pure  gases  and 
gaseous  mixtures  of  van  der  Waals  and  his  associates  at  the 
Universities  of  Amsterdam  and  Leiden. 

Figure  3.1  shows  the  p-v-T  relationship  of  the  liquid  and 
gas  phases  of  a  pure  fluid.  The  upper-right  projection  shows 
several  isotherms  on  a  p-v  diagram,  the  upper-left  projection 
is  a  p-T  plot  showing  several  isochores  and  the  bottom 
projection  shows  several  isobars  on  a  T-v  plot.  The  tongue- 
shaped  region  bounded  by  points  A,  D,  C,  E,  B  is  the 
coexistence  curve  (or  vapor  pressure  curve.)  This  curve  may  be 
considered  as  consisting  of  two  curves:  the  bubble-point  (the 
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Figure  3.1.  The  p-v-T  behavior  of  pure  fluid.  In  the  center  is 
sketched  the  surface  p  =  p(v,T) .  [From  Hirschfelder  Joseph  0., 
Curtiss  Charles  F.,  and  Bird  R.  Byron.  Molecular  Theory  of 
Gases  and  Liquids.  Copyright  ©  1954  by  John  Wiley  &  Sons,  Inc. 
Reprinted  by  permission  of  John  Wiley  &  Sons,  Inc.] 

point  of  initial  vaporization  when  the  pressure  of  the  liquid 
is  reduced)  curve,  ADC,  and  the  dew-point  (the  point  of 
initial  condensation  when  the  pressure  of  the  gas  is 
increased)  curve,  BEC.  These  two  curves  meet  at  the  critical 
point  C  which,  belonging  to  both  curves,  indicates  the 
identity  of  the  liquid  and  vapor  phases.  At  this  point, 
viewing  from  the  isothermal  perspective, 
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In  pure  fluid  systems  the  phenomenon  of  condensation  is 
associated  with  lowering  the  temperature  and  raising  the 
pressure  and  vice  versa  for  the  phenomenon  of  vaporization. 

In  a  fluid  mixture  system  we  have  besides  the  external 
parameters:  temperature  and  pressure,  the  internal  parameters 
which  identify  the  composition.  The  phase  behavior  of  the 
system  can  be  more  complicated  than  that  of  a  pure  fluid 
system  due  to  many  factors  such  as  the  possibility  of  various 
kinds  of  retrograde  phenomena,  the  occurrence  of  azeotropic 
mixtures  of  positive  or  negative  types  etc.  In  a  binary 
system,  for  example,  on  the  basis  of  an  analysis  of  the  phase 
diagrams  by  means  of  the  van  der  Waals  equation  of  state,  nine 
major  types  of  phase  diagram  may  be  obtained*42"44' .  The  complete 
phase  behavior  of  a  binary  system  is  represented  by  the  four 
dimensional  surface,  p  =  p(v,t,x)  where  x  is  the  mole  fraction 
of  one  of  the  two  components.  Consequently,  the  critical  point 
of  a  mixture  has  to  be  redefined  as  the  point  where  liquid  and 
vapor  become  identical,  subject  to  the  constraint  of  equality 
of  composition. 

Since  this  work  deals  with  a  binary  mixture  of  carbon 
dioxide  and  ethane,  attention  throughout  will  be  mainly 
focused  on  the  relevant  information  for  this  mixture.  Carbon 
dioxide  and  ethane  form  an  azeotropic  mixture  at  the 
appropriate  composition.  This  binary  system  deviates 
positively  from  the  ideal  solution,  which  obeys  Raoult's  law. 
In  this  case  a  critical  point  locus  follows  a  curve  between 
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Figure  3.2  Pressure-Temperature-Mole  fraction  relationship 
of  system  of  carbon  dioxide  and  ethane. 
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Figure  3.3.  Pressure-temperature  projection  of  the  system 
carbon  dioxide-ethane  in  the  critical  region.  Lines  Z  K 
and  ZbKb  are  vaP°r  pressure  curves  of  pure  carbon  dioxide 
and  of  pure  ethane  respectively. 
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the  critical  point  of  pure  carbon  dioxide  and  that  of  pure 
ethane.  Figures  3.2  and  3.3  illustrate  the  pressure- 
temperature-mole  fraction  diagram  and  critical  locus  of  this 
mixture,  respectively.  In  figure  3.3  at  point  B  the  maximum 
azeotrope  is  tangent  to  the  critical  curve  (KABKB)  ,  which  has 
a  minimum  in  temperature . 


Determinations  of  Critical  Point 

Conventionally,  there  are  three  principal  methods  of 
locating  the  critical  point.  These  are  described  briefly 
below.  (45> 
Isotherm  Approach 

Critical  temperature  may  be  estimated  to  within  a  few 
hundredths  of  a  kelvin  by  an  analysis  of  the  geometry  of  the 
isotherms'46'.  The  isotherms  in  the  immediate  neighborhood  of 
the  critical  temperature  are  measured  with  sufficient  accuracy 
to  allow  the  critical  temperature  to  be  determined  from  the 
inflection  points  of  the  isotherms.  The  minimum  values  of  the 
derivative  of  pressure  with  respect  to  volume  at  constant 
temperature,  \-^-\  found  graphically  are  plotted  as  a 
function  of  temperature  (or  pressure  or  molar  volume.)  The 
intersection  of  this  line  with  the  horizontal  axis  gives  the 
critical  parameters  (T.,  Pc,  vc.  )  Figures  3.4-3.6  illustrate 
this  approach. 
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Isobar  Approach 

This  method  is  sometimes  called  the  rectilinear  diameter 
method  of  Cailletet  and  Mathias.(48)  If  pL  and  p  are  the 
densities  of  liquid  and  of  saturated  vapor  in  equilibrium  with 
it  (so  called  orthobaric  densities),  their  mean  is  a  linear 
function  of  temperature. 


-j  (Pi  +  Pg)  "  P0  +  at 


(3.2) 


where  p0  =  mean  density  of  liquid  and  its  saturated  vapor  at 
0  °C 

t  ■  temperature  in  Celsius 

a  =  constant. 
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Figure  3.4.  Schematic  diagram  of  pressure  and  volume 
relationship  of  carbon  dioxide  with  several  isotherms  in 
a  broad  region. 
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Figure  3.6.  Schematic  diagram  of  graphical  determination  of 
critical  temperature  of  carbon  dioxide. 
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Thus,  if  p  and  p1  are  separately  plotted  against  the 
temperature  t,  the  locus  of  the  points  bisecting  the  joins  of 
corresponding  values  of  p1  and  p  is  a  straight  line  (see 
figure  3.7.)  The  point  where  this  straight  line  cuts  the 
coexistence  curve  is  the  critical  point.  In  many  cases  this 
empirical  law  of  rectilinear  diameters  holds  very  well,  but 
sometimes  over  a  large  range  of  temperature  the  actual 
behavior  shows  a  slight  curvature.  Furthermore,  very  near  the 
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Cailletet  and  Mathias. 
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critical  point  (within  several  mK  of  Tc)  real  systems  show  a 
deviation  from  linearity. 

Isochore  Approach 

There  are  two  slightly  different  ways  for  taking  this 
approach.  The  first  one  is  to  study  the  discontinuity  of  the 
isochore  curve.'"1  This  method  begins  with  loading  a  bomb  of 
constant  volume  with  a  series  of  known  weights  of  the 
substance  and  then  studying  the  behavior  of  the  pressure  as 
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Figure  3.8.   Schematic  diagram  of  isochoric  method  of 
critical  point  determination. 
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the  temperature  changes.  The  following  possible  kinds  of 
isochoric  curves  could  result  if  the  charging  density  p,is 
greater  than  the  critical  density  pc,  the  isochore  curve  shows 
an  upward  bend  as  in  C,  or  if  p  is  less  than  pc  the  curve 
shows  downward  bend  as  in  A.  However,  if  p  is  equal  to  p  the 
isochore  exhibits  no  break  as  in  B.  This  method  is  less 
demanding  in  terms  of  skillful  technique  and  special 
apparatus,  but  the  accuracy  is  poor.  Figure  3.8  shows  data 
illustrating  all  three  possibilities  described  above. 

The  second  established  isochoric  method  is  to  study  the 
appearance  and  disappearance  of  the  meniscus  between  the 
phases  near  the  critical  temperature. <19)  The  procedure  is  the 
same  as  above  except  that,  instead  of  measuring  pressure  as 
temperature  changes,  the  behavior  of  the  meniscus  separating 
the  liquid  and  gaseous  phases  is  observed.  There  are  also 
three  possible  phenomenon: 

If  p  <  pc  the  meniscus  falls  until  the  entire  container 
is  filled  with  gas. 

If  p  >  pc  the  meniscus  rises  until  the  entire  container 
is  filled  with  liquid. 

If  p  =  pc  the  meniscus  approximately  located  at  a  point 
halfway  up  the  container  will  flatten,  then 
become  very  faint  and  finally  disappear. 
The  temperature  and  density  at  which  the  meniscus  disappears 
are  taken  to  be  the  critical  parameters. 
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Speed  of  Sound  as  a  Probe  for  the  Critical  Point 

By  conventional  methods  the  critical  point  can  be  located 
indirectly  as  the  zero  point  of  the  first  derivative  of 
pressure  with  respect  to  volume  or  of  volume  with  respect  to 
temperature  on  an  isotherm  or  isobar,  respectively,  or  by 
analyzing  the  slope  of  lines  plotted  as  the  first  derivative 
of  pressure  with  respect  to  temperature  versus  temperature  for 
isochores.  The  technique  involving  use  of  the  principle  of 
appearance  and  disappearance  of  a  meniscus  is  to  some  degree 
subjective  due  to  the  heuristic  effect.  A  better  alternative 
would  be  to  use  some  thermodynamic  property  which  shows  a 
significant  change  at  the  critical  point.  In  1962  Bagatski  and 
his  coworkers  <50'51'  experimentally  observed  the  asymptotic 
behavior  of  heat  capacity  at  constant  volume  of  argon  in  the 
immediate  vicinity  of  the  critical  point.  As  expected  C  falls 
off  rapidly  and  it  appeared  that  this  could  be  used  to  detect 
the  critical  point.  Measuring  heat  capacity,  however,  is  not 
an  easy  task.  On  the  other  hand,  the  speed  of  sound  similarly 
shows  an  abrupt  change  toward  zero  at  the  critical  point  since 
it  is  inversely  proportional  to  Cv  which  diverges  weakly  and 
is  directly  proportional  to  C  which  diverges  strongly  and  to 
the  first  derivative  of  pressure  with  respect  to  molar  volume 
at  constant  temperature,  which  also  vanishes  at  the  critical 
state.  The  appropriate  equation  is: 
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where  C0  is  the  sonic  speed  in  the  limit  of  zero  frequency. 
Because  the  sonic  speed  is  directly  proportional  to  the 
frequencies  of  normal  mode  resonance  vibrations  in  a  filled 
cavity  resonator,  one  can  expect  to  detect  the  approach  to  the 
critical  state  by  observing  changes  in  resonance  frequency. 
Spherical  cavities  yield  especially  sharp  resonances  and  are 
therefore  well  suited  for  this  purpose.  In  the  first  phase  of 
this  work  the  isochoric  approach  was  employed.  However,  unlike 
the  conventional  method  of  using  the  straightness  of  the 
pressure-temperature  plot  or  using  the  appearance  and 
disappearance  of  the  meniscus,  the  speed  of  sound  or  a 
resonance  frequency  was  utilized  as  a  probe  of  a  critical 
point. 


CHAPTER  4 
EXPERIMENT 


Instruments 


The  development  of  a  spherical  acoustic  resonator  as  a 
tool  to  the  study  speed  of  sound  in  gases  such  as  argon, 
butane,  isobutane  etc.  which,  in  turn,  leads  to  values  of  some 
thermodynamic  and  equation  of  state  properties  has  been 
successful  in  the  past  by  workers  in  this  research  group.  (20"26) 
There  are,  however,  some  limitations  on  the  previous  work  due 
mainly  to  the  construction  of  the  spherical  resonator.  It  was 
made  from  aluminum  alloy  and  designed  to  operate  near 
atmospheric  pressure.  These  constraints  limit  its  applications 
substantially.  For  example,  if  one  wishes  to  study  natural  gas 
mixtures,  one  would  likely  encounter  the  presence  of  some 
corrosive,  acidic  or  basic  gases  such  as  hydrogen  disulfide, 
carbon  dioxide  etc.,  and  these,  especially  in  the  presence  of 
water,  attack  aluminum  to  some  extent.  In  addition, 
experiments  may  call  for  the  measurements  at  much  higher 
pressures  than  atmospheric. 

The  present  research  was  undertaken  to  extend  the 
applicability  of  the  acoustic  resonance  technique  to  the  study 
of  natural  gas  mixtures  at  high  pressures.  To  fulfil  this  goal 
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a  new  high  performance  stainless  steel  resonator  assembly  was 
designed  and  built,  and  new  data  acquisition  software  was 
developed. (52)  The  new  apparatus  was  used  to  measure  the  sonic 
speed  in  some  carefully  blended  gas  mixtures.  Because  the 
mixtures  contained  small  concentrations  of  condensibles  such 
as  hexane  for  example,  it  was  possible  to  cool  them  below 
their  dew  points  and  it  was  observed  that  the  sonic  speed 
faithfully   revealed  the  precise   location   of   the   phase 
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Figure  4.1.  Observed  sonic  speed  versus  temperature  for  a 
North  Sea  natural  gas  mixture. (Reference  53) 
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change'53'.  This  behavior  was  soon  recognized  as  a  valuable 
tool  for  phase  equilibrium  studies.  Figure  4.1  shows  the 
typical  behavior  in  the  variation  of  sonic  speed  with 
temperature  for  such  a  mixture.  The  dew  point  is  clearly 
evident.  Following  this  discovery  it  was  decided  to  evaluate 
the  use  of  the  acoustic  resonator  to  locate  the  most 
interesting  feature  of  the  phase  diagram,  namely  the  critical 
point.  The  work  reported  in  this  dissertation  deals  primarily 
with  this  difficult  and  challenging  task. 
Spherical  resonator 

The  resonator  cavity  is  the  heart  of  the  apparatus.  To 
permit  its  use  at  high  pressure  and  in  the  environment  of 
corrosive  gases,  the  spherical  resonator  was  built  from 
stainless  steel.  Its  side  view  is  shown  in  figure  4.2. 

The  resonator  cavity  was  fashioned  from  a  welded 
spherical,  type  304  stainless  steel  tank  approximately  0.203 
m.(8  in.)  in  diameter  with  a  0.24  cm. (3/32  inch.)  thick  wall. 
The  tank  is  not  perfectly  spherical.  This  imperfection, 
however,  does  not  significantly  degrade  the  data  integrity  as 
long  as  appropriate  measurement  techniques  are  used  as 
discussed  in  the  pervious  section  and  proven  to  be  valid  by 
the  previous  investigators.20'21  The  tank  was  obtained  from  the 
manufacturer  (Pollution  Measurement  Corp. , Chicago,  IL)  with 
only  a  single  0.6  cm.(%  in.)  hole  located  at  one  pole  (6=0 
in  spherical  polar  coordinates).  A  matching  hole  was  drilled 
through  the  opposite  pole  (0  =  n)    and  two  0.953  cm  (%  in.) 
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Figure  4.2.  A  side  view  of  the  spherical  resonator  equippec 
with  the  transducers. 
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Figure  4.3.  The  spherical  resonator. 
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holes  drilled  through  the  wall  at  0  =  tt/4  with  <f>  =  0  and  n, 
respectively.  No  holes  were  located  on  the  equatorial  weld 
bead  (8  =  n/2)  .  Figure  4.3  shows  the  overall  view  of  the 
spherical  resonator. 
Transducers 

The  active  element  of  the  acoustic  transducer  is  a  0.750 
in.  diameter  piezoelectric  ceramic  bimorph  (Vernitron, 
Bedford,  Ohio)  made  of  PZT-5A  material  (lead  zirconate-lead 
titanate.)  This  ceramic  was  found  to  produce  a  signal  of  high 
volume  when  driven  with  a  sinusoidal  waveform  of  10  volts 
peak-to-peak.  It  also  has  a  high  maximum  operating  temperature 
of  250  °C.(54)  A  transducer  assembly  was  designed  and 
fabricated.  It  is  shown  in  figure  4.4.  The  electrical 
feedthrough  made  from  brass  was  held  in  a  Delrin  insert  to 
prevent  a  short  circuit.  A  cross-sectional  view  of  the 
transducer  mounted  on  the  resonator  is  also  shown  in  figure 
4.2.  The  transducer  assemblies  were  mated  to  the  resonator 
through  short  %  in.  diameter  tubes. 

Two  identical  transducers  located  at  right  angles  to  each 
other  were  used.  This  orientation  has  proved  to  yield  a  better 
resolution  of  the  radial  modes  of  vibration  due  to  diminishing 
in  intensity  of  some  non-radial  modes,  p  modes  for  example. 
One  transducer  functions  as  a  speaker.  It  is  an  input 
transducer  which  mechanically  deforms  when  a  voltage  is 
applied.  The  other  transducer  functions  as  a  microphone.  It  is 
an  output  transducer  which  operates  on  a  reverse  mechanism  of 
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Figure  4.4.  The  transducer  (cross-sectional  view). 
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the  input  transducer:  producing  an  electric  signal  when 
mechanically  deformed  by  pressure  waves. 
Mixing  control  system 

A  circulating  pump  is  needed  to  promote  mixing  of  system 
components  and  (or)  phases  and  to  assist  with  thermal 
equilibration  by  eliminating  stratification.  It  was  made  from 
stainless  steel  and  consists  of  two  parts:  the  liquid 
collector  and  the  pump.  Figure  4.5  shows  the  design  of  the 
unit.  It  consists  of  a  vertically  mounted  stainless  steel 
cylinder  3.18  cm  (1.25  in.)  in  diameter  with  a  freely  sliding 
internal  piston  2.10  cm  (0.851  in.)  in  diameter.  The  piston 
has  a  built  in  reed  valve  which  is  normally  open  under  gravity 
but  is  forced  closed  by  viscous  drag  as  the  piston  descends 
causing  the  fluid  to  flow  from  top  to  bottom  through  the  pump 
on  descending  motion  of  the  piston.  Piston  motion  is  driven  by 
interaction  between  an  iron  slug  attached  to  the  piston  and  an 
external  magnet  driven  in  turn  by  linkage  to  a  motorized  bell 
crank. 
Apparatus  assembly 

Tubular  fittings  to  the  spherical  resonator  to  other 
parts  of  the  system  were  provided  with  0.6  cm  (h  in.)  diameter 
pilot  tubulations  which  slip  into  the  polar  holes  and  when 
pressed  against  the  tank  were  sealed  by  elastomeric  O-rings 
confined  to  glands  comprised  of  circular  grooves  in  the 
fittings  around  the  pilot  tubes  and  the  ring  like  areas  of  the 
adjacent  tank  wall  surrounding  the  holes.  These  fittings  were 
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Figure   4.5.    The  mixing  control   unit. 


41 
brazed  to  stainless  steel  tubes  used  for  charging  the 
resonator  and  circulating  its  contents  through  an  external 
loop  containing  the  magnetically  driven  pump.  The  two 
identical  PZT  bimorph  piezoelectric  transducer  assemblies  were 
similarly  mounted  to  the  tank  at  the  0.953  cm  (%  in.)  hole 
positions . 

The  assembly  support  frame  consists  of  two  parallel  0.6 
cm  (h  in.)  thick  stainless-steel  plates  between  which  an  array 
of  2.54  cm  (1  in.)  thick  stainless-steel  blocks  were  attached 
with  bolts  to  support  the  tube  fittings  and  transducer 
assemblies  and  provide  a  means  of  applying  clamping  forces  to 
the  O-ring  seals.  The  clamping  force  acting  along  the  polar 
axis  is  applied  by  a  hallow  jam  screw  which  slips  over  the 
charging  tube  and  pushes  against  the  brazed  O-ring  gland.  The 
block  which  supports  the  lower  tube  fitting  is  rigidly  fixed 
and  provides  a  definite  reproducible  location  for  the  tank 
body.  The  transducer  assemblies  were  clamped  against  the 
resonator  by  set  screws  threaded  through  two  of  the  blocks 
bolted  to  the  mounting  frame.  Set  screws  in  two  additional 
blocks  located  opposite  the  transducer  mounting  holes  apply 
counter  forces  against  the  tank  to  stabilize  the  assembly 
further. 

The  circulating  pump  is  connected  to  the  upper  and  lower 
tube  fittings  by  demountable  compression  fittings  (Swagelok.) 
The  normal  pumping  rate  of  4  liters  per  minute  proved  adequate 
for  these  purposes.  A  0.6  cm  (%  in.)  tube  connected  to  the 
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fluid  circulation  loop  by  a  tee  leads  to  a  diaphragm-type 
pressure  gauge  and  a  system  shut-off  valve. 

The  resonator  assembly  and  circulating  pump  were 
suspended  from  a  horizontal  plate  which  forms  the  cover  of  a 
stirred  liquid  bath  in  which  the  system  is  mounted.  The  bath 
was  housed  in  a  well-insulated  container  and  connected  by 
insulated  tubing  to  a  computer-controlled  heat  exchanger 
capable  of  operating  from  233  to  373  K.  Bath  temperature  was 
measured  with  a  four-wire  platinum  resistance  thermometer 
(RTD.)  Figure  4.6  shows  the  instrumental  setup. 
Electronic  hardware 

The  automated  systems  employed  in  this  work  are  shown  in 
figures  4.7,  4.8,  and  4.9.  Figure  4.7  shows  the  setup  used  in 
the  experiments  on  pure  carbon  dioxide  and  on  a  mixture  of 
carbon  dioxide  and  ethane.  Figure  4.8  shows  the  setup  for  the 
experiments  on  pure  ethane.  In  this  latter  setup  we  have 
explored  the  capability  of  using  a  fast-fourier  transform 
(FFT)  technique.  Unfortunately,  the  data  collected  by  this 
technique  was  unreliable  due  to  leaks  in  the  system  found 
after  the  experiment  was  finished.  While  the  leaks  were  being 
fixed,  the  instruments  were  transferred  to  another  project. 
Consequently,  the  FFT  exploration  was  discontinued  for  ethane 
and  the  setup  shown  in  figure  4.9  was  built  to  finish  this 
project. 
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Even  though  there  were  different  setups,  the  basic 
principles  are  the  same  and  can  be  divided  into  two  parts :  the 
input  signal  generation  part  and  the  output  signal  analysis 
part.  These  are  shown  by  dashed-line  boxes  in  figures  4.7  and 
4.9.  Commands  and  data  were  transferred  among  the  instruments 
and  computer  over  the  General  Purpose  Interface  Bus  (GPIB) 
cables.  Each  component  on  the  GPIB  was  equipped  with  an  IEEE- 
488  standard  interface.  The  computer  was  always  a  controller- 
in-charge  sending  commands  and  acquiring  data  according  to  the 
written  computer  program.  Each  instrument  was  a  listener 
and/or  a  talker  depending  on  a  currently  executed  line  of  the 
program.  Table  4.1  summarizes  electronic  hardware  components 
used  in  this  work  and  their  specifications. 

Gases 

Carbon  dioxide  and  ethane  gases  used  were  high  purity 
grade.  No  further  analysis  was  performed  on  carbon  dioxide. 
Table  4.2  summarizes  the  pure  gas  qualities. 

For  the  experiments  on  a  mixture  of  C02-C2H  quantitative 
analyses  were  performed  since  the  results  of  speed  of  sound 
measurements  depend  on  the  precise  composition.  Both 
gravimetric  and  gas  chromatography  analyses  were  used  to 
fulfil  this  task.  Details  of  these  analyses  are  presented  in 
appendix  B. 
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Table  4.2.  Gas  Specifications. 

Gas 

Manufacturer 

Grade 

Purity 

Carbon  dioxide 

Scott 

specialty 

gases 

Research 
Grade 

99.99  Mole  % 

Ethane 

Scott 

Specialty 

Gases 

Research 
Grade 

99.9  % 

Argon 

Matheson  gas 
Products 

Research 
Grade 

99.9995  Volume  % 

Table  4.3.  Charging  Pressure  of  gas  systems, 


Gas  system 

Starting 
temperature  (°C) 

Charging  pressure 
(psia) 

Carbon  dioxide 

36.0 

1200 

C02   +  C2H6 

30.0 

630  for  C2H6  and 
make  up  to  1300 
with  CO, 

cA 

37.0 

808 

Resonance  Frequency  Measurements 


Measurements  were  made  on  isochore  samples  introduced 
into  the  spherical  resonator  which  was  initially  rinsed  with 
argon  gas  and  subsequently  evacuated  overnight  at  high 
temperature  (40-50  °C)  with  a  mechanical  pump  connected  to  a 
liquid  nitrogen  trap.  The  apparatus  was  then  brought  to  a 
temperature  above  the  known  literature  critical  temperature. 
After  that  gas  was  introduced  slowly  into  the  resonator  to 
bring  the  pressure  to  the  calculated  value  (see  table  4.3). 
The  apparatus  was  then  sealed  off.  The  charging  pressure  value 
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was  calculated  using  the  AGA8  program.  Details  of  the 
calculation  are  given  in  appendix  A.  In  the  case  of  the 
mixture  of  C02  and  C2Hg  ethane  was  first  introduced.  Pressure 
was  monitored  roughly  by  a  Bourdon  pressure  gauge  and  more 
precisely  by  a  digital  pressure  gauge  previously  calibrated  to 
yield  an  absolute  pressure  reading.  The  magnetic  circulating 
pump  was  turned  on  to  assure  adequate  mixing.  Once  the  system 
was  at  stable  conditions  in  temperature  and  pressure  (which 
took  approximately  a  few  hours),  searching  for  resonance 
frequencies  began. 

Fluids  in  the  resonator  were  stimulated  acoustically  by 
driving  one  of  the  bimorph  transducers  with  a  40  volt  peak-to- 
peak  sine-wave  signal  generated  under  microprocessor  control 
(PC-AT)  of  a  programmable  signal  synthesizer,  HP3325B  in  the 
case  of  pure  carbon  dioxide  and  of  the  of  C02  and  C2H6  mixture 
or  of  a  wave  function  generator  (WAVETEX)  controlled  through 
a  programmable  digital-to-analog  converter  (DAC)  in  the  case 
of  pure  ethane.  The  fluid  response  was  determined  at  the  other 
transducer  which  feeds  its  signal  to  a  lock-in  amplifier.  The 
response  of  the  speaker  transducer  in  terms  of  voltage  output 
signal  was  monitored  closely  to  changes  in  frequency  of  the 
microphone  transducer.  At  resonance,  the  output  amplitude 
increased  dramatically  to  maximum  enabling  a  precise  frequency 
measurement.  Using  an  oscilloscope  as  a  visual  aid  to  locate 
the  vicinity  of  this  maximum  signal  and  a  digital  voltmeter  to 
find  the  precise  position  of  the  signal  peak,  resonance 
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frequency  could  be  found  with  high  accuracy.  Usually,  the 
first  radial  mode  of  vibration  was  selected  as  the  tracked 
peak  throughout  this  work. 

The  approximate  location  of  the  desired  peak,  first 
radial  mode  for  example,  was  initially  calculated  from  the 
AGA8  equation  using  pressure,  temperature,  and  composition  as 
input.  Then  the  actual  location  of  the  peak  was  experimentally 
searched  using  a  program  called  SCAN.  Scanning  was  performed 
in  a  reasonable  frequency  range  around  the  center  frequency 
identified  as  the  output  value  of  the  AGA8  calculation.  Once 
the  peak(s)  was  found,  location  of  another  mode  was  calculated 
as  described  in  the  theoretical  background  section.  The 
purpose  of  this  routine  was  to  increase  the  degree  of 
confidence  that  each  experimentally  found  peak  was  indeed 
correctly  identified.  Similarly,  if  more  than  one  peak 
appeared  together  in  the  range  while  scanning,  each  was 
checked  against  the  other  radial  modes  to  sort  out  the  best 
candidate  for  labeling  as  a  resonance  mode  based  on  the 
hypothesis  that  if  it  is  a  genuine  system  peak,  its  frequency 
should  reveal  the  correct  position  of  the  other  peaks.  Figure 
4.11  shows  an  example  of  a  series  of  scans  performed 
chronologically.  The  figure  in  the  first  row  is  an  acoustic 
spectrum  scanned  in  a  wide  frequency  range.  The  figure  in  the 
second  row  is  the  acoustic  spectrum  scanned  in  a  narrower 
frequency  than  the  one  in  the  previous  figure.  In  this 
incident  there  were  three  peaks  near  the  predicted  frequency 
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of  the  first  radial  mode  peak  (peak  1  in  figure  4.11).  Each  of 
these  was  investigated  closely  as  shown  in  the  subsequent  four 
rows.  Each  figure  in  the  third  row  is  a  spectrum  of  a  single 
peak  of  row  two  scanned  in  a  much  narrower  frequency  range. 
Frequencies  found  for  each  peak  in  this  row  were  used  to 
calculate  the  second,  third,  and  forth  radial  frequencies. 
Peaks  corresponding  to  these  calculated  frequencies  were 
searched  as  shown  the  figures  in  the  forth,  fifth,  and  sixth 
row.  Apparently,  the  third  peak  (peak  3  in  figure  4.10)  is  not 
the  correct  peak  as  it  failed  to  predict  the  location  of 
higher  radial  modes.  Peaks  1  and  2  give  similar  predictions 
since  they  lie  very  close  together;  about  5  Hz  apart. 
Predictions  of  overtones  based  on  peak  1  were  in  better 
agreement  with  experiment  than  from  peak  2  therefore  peak  1  is 
correctly  interpreted  as  the  first  radial  mode  peak.  Note  that 
the  SCAN  program  was  used  for  pure  carbon  dioxide  and  for  the 
C02-C2H6  mixture  experiments.  For  pure  C2Hfi  experiment  this 
process  was  performed  manually  via  a  voltage-control-generator 
(VCG)  box  since  the  instrumental  setup  was  altered. 

Having  established  the  location  of  the  peak  to  be 
tracked,  data  collection  was  begun.  For  pure  CO  and  for  the 
mixture  of  C02-C2Hfi  experiments  programs  called  VIKING  and 
VIKAN  were  used.  In  this  version  only  data  corresponding  to 
the  maximum  amplitude  were  recorded  for  each  sweep  frequency 
range.  For  the  pure  C2H6  experiment,  however,  data  of  the 
entire  sweep  frequency  range  were  recorded  (via  the  program 
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Figure  4.10.  Series  of  scanning  routine, 
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called  MAX)  since  it  was  believed  the  more  detailed 
information  would  be  instructive.  This  change  was  possible 
because  of  advances  in  the  available  data  acquisition 
hardware.  In  both  cases  the  programs  locked  onto  the  tracked 
peak  for  the  entire  process  as  bath  temperature  was  lowered 
(or  raised)  at  a  rate  of  a  few  millikelvin  per  minute  over  a 
period  of  several  hours.  Pressure  and  temperature  measurements 
were  simultaneously  collected  along  with  frequency  and 
amplitude.  Upon  completion  of  the  experiment  the  temperature 
of  the  system  was  brought  back  to  the  starting  value.  A  small 
amount  of  gas  was  removed  to  change  the  density  of  the  system, 
and  then  another  data  collection  was  begun.  Typically,  the 
density  of  the  system  was  initially  at  a  supercritical  value 
(p  >  pc)  and  was  finally  at  a  subcritical  value  (p  <  p  )  .  The 
usual  number  of  runs  for  a  system  was  about  fifty. 

Note  also  that  the  resonance  frequency  may  also  be 
determined  accurately  by  observation  of  the  phase  change  of 
the  output  voltage  relative  to  the  excitation.  This  approach 
is  necessary  for  the  most  accurate  determination  of  absolute 
sonic  speeds,  but  great  care  must  be  taken  since  the  phase 
meter  is  extremely  sensitive  and  responds  to  peaks  at  even  the 
noise  level.  All  programs  were  written  in  FORTRAN  77.  They  are 
shown  in  appendix  E.  Figure  4.11  shows  typical  output  response 
of  these  two  approaches. 

Once  data  collection  was  complete  the  corresponding  speed 
of  sound  were  calculated  with  a  computer  program  called 
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GETSPG.  This  program  utilizes  the  principle  of  relative 
measurement  mentioned  earlier  in  the  theory  section.  Also  the 
theoretical  speeds  of  sound  were  calculated  from  the  AGA8 
equation  through  a  computer  program  call  TOTDAT  for  comparison 
purposes.  Finally  the  results  were  manipulated  and  graphically 
displayed  using  popular  software  including:  LOTUS-123, 
QUATTRO,  GRAFTOOL,  and  GRAPHER.  Figure  4.12  shows  a  schematic 
diagram  of  the  entire  process  of  data  collection. 
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CHAPTER  5 
RESULTS  AND  DISCUSSIONS 


Carbon  Dioxide 


As  mentioned  in  the  introduction,  due  to  its  well- 
established  values,  carbon  dioxide  was  chosen  to  be  the  first 
candidate  used  to  test  this  new  acoustic  resonance  technique 
for  locating  critical  point  parameters.  Learning  how  best  to 
do  this  required  considerable  effort,  often  involving  trial 
and  error  procedures.  Many  adjustments  had  to  be  made  to  both 
instruments  and  programs  to  identify  the  optimum  conditions 
for  performing  this  experiment. 

Figure  5 . 1  shows  the  typical  relationship  between  the 
first  radial  mode  resonance  frequency  (hereafter  abbreviated 
resonance  frequency)  and  the  temperature  of  carbon  dioxide  in 
the  resonator  charged  near  its  critical  density.  It  shows  the 
results  of  four  runs  of  this  same  isochore.  For  each  run  the 
curve  may  be  roughly  divided  into  two  parts  for  the  sake  of 
discussion.  The  first  one  is  a  curve  where  the  system 
temperature  is  greater  than  or  equal  to  the  critical 
temperature,  Tc.  In  this  region  the  speed  of  sound  is  directly 
proportional  to  the  system  temperature,  i.e.  as  temperature 
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Figure  5.2.  Relationship  between  first  radial  mode  resonance 
frequency  and  temperature  of  carbon  dioxide  for  an  isochore 
near  the  critical  density. 


Figure  5.3.  Relationship  between  the  speed  of  sound  and 
temperature  of  carbon  dioxide  for  an  isochore  near  the 
critical  density. 
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decreases,  resonance  frequency  and,  hence  speed  of  sound, 
decreases.  This  represents  a  single  phase  region  of  the  fluid. 
The  second  part  is  a  curve  where  the  system  temperature  is 
less  than  the  critical  temperature.  Unlike  the  first  region, 
the  speed  of  sound  is  inversely  proportional  to  the  system 
temperature,  i.e.  as  temperature  decreases,  resonance 
frequency  increases.  This  represents  a  two-phase  region  where 
both  gas  and  liquid  are  present  in  equilibrium.  Each  of  the 
experimental  plots  has  a  broken  tip  which  occurs  in  the 
critical  region  indicating  loss  of  signal  and  hence  loss  of 
precise  track  of  the  resonance  peak  in  this  vicinity.  This  is 
understandable  in  view  of  the  extremely  fast  dynamic  behavior 
of  the  resonance  frequency  in  the  critical  vicinity  and  the 
simultaneous  dramatic  reduction  in  signal  amplitude.  The  best 
experimental  result  was  obtained  during  a  run  in  which  the 
temperature  of  the  well  insulated  system  was  initially  set  at 
a  value  above  the  estimated  critical  temperature,  and  the 
system  was  then  allowed  to  cool  down  naturally  to  room 
temperature  (around  296  K.)  The  cooling  time  was  approximately 
17  hours  during  which  approximately  1000  data  points  were 
collected.  Figures  5.2  and  5.3  represent  results  of  this 
particular  run  showing  the  relationships  between  system 
temperature  and  resonance  frequency  and  sonic  speed 
respectively.  Results  of  this  run  clearly  show  the  extremely 
sensitive  dependence  of  sonic  speed  on  temperature  in  the 
critical  region.   Figure  5.4  shows  a  plot  of  the  first 
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Figure  5.5.  Enlarged  temperature  scale  of  figure  5.4  showing 
the  critical  temperature  to  be  at  304.215  K. 
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derivative  of  resonance  frequency  with  respect  to  temperature 
versus  system  temperature  for  the  data  of  figure  5.3.  It  gives 
a  precise  value  of  the  critical  temperature  as  seen  in  figure 
5.5  which  is  a  portion  of  figure  5.4  with  an  expanded 
temperature  scale.  The  critical  temperature  obtained  from  this 
plot  is  304.215±0.007  kelvin.  Figure  5.6  displays  the 
relationship  between  the  system  pressure  and  temperature  for 
this  same  run.  It  shows  a  nearly  straight  line  with  no 
evidence  of  phase  separation. 

Experiments  were  also  run  with  subcritical  density 
charges.  These  show  the  phase  boundary  location  (dew  points) 
as  changes  in  slope  of  resonance  frequency  versus  system 
temperature.  Table  5.1  gives  a  comparison  of  results  found  by 
the  present  method  with  available  literature  values.  The 
present  results  are  in  good  agreement  with  the  other 
experimental  values  and  with  predictions  based  on  an  equation 
of  state(55)  and  a  correlation  model  developed  at  the  National 
Institute  of  Standards  and  Technology  (NIST,  formerly  NBS)<56). 
The  acoustic  data  and  NIST  predictions  are  shown  graphically 
in  figure  5.7.  Table  5.2  gives  a  chronological  collection  of 
experimental  values  of  critical  point  parameters  (T  and  P  ) 
reported  by  several  authors.  The  value  obtained  by  this  sonic 
method  is  only  0.0016  percent  different  from  the  best  value 
obtained  in  the  NIST  laboratories  (Morrison  and  Kincaid,  1984. 
Reference  96) . 
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Figure   5.7.    Experimental    vapor    pressure    curve    compared   with 
the  NIST  model. 
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Table  5.2.  Values  of  critical  point  parameters  from  the 
literature. 


Temperature 
K 

Pressure 

Density 

Author (s)* 

atm. 

psia. 

Kg/m3 

304.07 

73.0 

1072.77 

- 

Andrews  (58) 

304.15 

- 

- 

- 

Hautefeuille 
Cailletet  (59) 

305.05 

77.0 

1131.55 

- 

Dewar  (60) 

304.50 

72.9 

1071.30 

- 

Amagat  (61) 

304.55 

- 

- 

- 

Chappuis  (62) 

304.85 

- 

- 

- 

Villard  (63) 

304.15 

- 

- 

- 

Verschaffelt  (64) 

304.25 

73.26 

1076.59 

- 

Kuenen  (65) 

304.55 

- 

- 

- 

De  Heen  (66) 

304.10 

- 

- 

- 

Von  Wesendonck 
(67) 

304.13 

72.93 

1071.74 

- 

Keesom  (68) 

304.27 

- 

- 

- 

Brinkmann  (69) 

304.135 

- 

- 

- 

Onnes  &  Fabius 
(70) 

304.41 

- 

- 

- 

Bradley,  Brown,  & 
Hale  (71) 

304.25 

73.00 

1072.77 

- 

Dorsman  (72 ) 

304.15 

72.85 

1070.56 

- 

Cardoso  &  Bell 
(73) 

304.12 

- 

- 

- 

Hein  (74) 

304.35 

- 

- 

— 

Dieterici  C.  (75) 

304.10 

- 

- 

- 

Pickering  S.F. 
(76) 

304.25 

72.95 

1072.03 

- 

Meyers  &  Van 
Dusen  (77) 

304.11±0.01 

— 

™ 

— 

Kennedy  H.T. , 
Cyril  H.,  and 
Meyer  (78) 

Table  5.2.  (Continued) 


69 


Temperature 
K 

Pressure 

Density 

Author (s)* 

atm. 

psia. 

Kg/m3 

304.11 

- 

- 

- 

Kennedy  (79) 

304.19 

— 

— 

— 

Michael  A. 
Blaisse  S. 
Michael  C.  (80) 

304.19 

- 

- 

- 

Lorentzen  and  Han 
Ludvig  (81) 

304.19 

72.80 

1069.83 

- 

Tielsch  H.  (82) 

304.16 

75.20 

1105.10 

- 

Ernst  S.  and 
Thomas  W.  (83) 

304.21 

72.87 

1070.86 

- 

Wentorf  R.H.  (84) 

304.16 
±0.03 

- 

- 

- 

Ambrose  D.  (85) 

304.15 
±0.01 

75.282 

1106.30 

- 

Ernst  S.  and 
Traube  K.  (86) 

304.2 

72.83 

1070.34 

468 

Matthews  J.F. 
(87) 

304.150 
±0.004 

72.79 

1069.68 

467.8 

Moldover  (88) 

304.19 

72.86 

1070.78 

468 

Altunin  V.  (89) 

304.13 

- 

- 

- 

Krynicki  (90) 

304.16 

- 

- 

- 

Lesnevskaya  (91) 

304.13 

72.79 

1069.62 

467 

Sengers  J.V.  (92) 

304.20 

- 

- 

- 

Morrison  G.  (93) 

304.13 

- 

- 

467.4 

Adamov  (94) 

304.13 

72.79 

1069.66 

467.8 

Shelomentsev  (95) 

304.206 
±0.001 

72.87 
±0.01 

1070.96 
±0.01 

468.248 
±0.008 

Morrison  Kincaid 
(96) 

304.215 
±0.007 

73.28 
±3.00 

1076.88 
±3.00 

- 

This  work 

Numbers  in  parenthesis  refer  to  references  in  bibliography. 
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Encouraged  by  the  apparent  success  of  this  approach  to 
the  location  of  phase  boundaries,  including  critical  points, 
for  a  pure  substance,  we  opted  to  proceed  to  the  next  stage  of 
the  study  and  perform  similar  tests  an  a  simple  fluid  mixture. 
The  properties  of  mixtures  are  not  so  well  documented  as  those 
of  carbon  dioxide,  and  there  are  few  previously  well 
characterized  binary  systems  to  choose  from.  Carbon  dioxide- 
ethane  is  one  of  these,  and  a  mixture  of  these  two  was 
selected  for  this  study.  The  results  obtained  are  presented  in 
the  next  section. 

A  Mixture  of  Carbon  Dioxide  and  Ethane 

Since  an  acoustic  resonance  technique  has  been 
successfully  used  to  detect  phase  changes  along  the 
coexistence  curve  of  a  single  component  system,  it  is  a 
naturally  logical  forward  step  to  investigate  its  suitability 
for  use  with  more  complex  systems.  What  is  needed  are  specific 
mixtures  with  well  characterized  phase  behavior.  A  binary 
mixture  of  74.25  mole  percent  carbon  dioxide  and  ethane  was 
selected  for  this  purpose.  The  carbon  dioxide-ethane  system 
has  been  studied  by  several  researchers  beginning  in  1897.  The 
first  critical  measurements  were  carried  out  by  Kuenen  in 
1897(97)  and  1902(98)  followed  by  Khazanova  and  Lesnevsicaya  in 
1967(99),  Fredenslund  and  Mollerup  in  1974<100»  and  Morrison  and 
Kincaid  in  1984(96).  Phase  boundary  studies  have  been  performed 
by  Clark  and  Din  in  1953(101),  Jensen  and  Kurata  in  1971(102), 


71 
Gugnoni  and  Eldridge  in  1974(103),  Nagahama  et  al.  in  1974(104), 
Dabalos  in  1974(105),  and  Stead  and  Williams  in  1980(106).  In  1978 
Moldover  and  Gallagher*107'  presented  a  correlation  for  this 
mixture  in  an  analogy  with  pure  fluids.  The  specific  mixture 
selected  for  this  study  is  one  which  was  studied  by  Khazanova 
and  Lesnevsicaya'99'.  Figures  5.8,  5.9,  and  5.10  show  the 
behavior  of  the  resonance  f reguency  versus  temperature  of  this 
mixture  at  three  different  gas  densities:  supercritical 
density  ( p  >  pc)  ,  subcritical  density  (p  <  p  ),  and  near 
critical  density  (p  aB  ~  pc) .  Note  the  sharper  turning  point  in 
the  curve  around  the  critical  region  in  figure  5.10  compared 
to  the  turning  points  of  the  curves  of  figure  5.8  and  5.9.  The 
former  is  a  distinctive  characteristic  of  a  system  at  its 
critical  conditions  (Tc,  Pc,  and  pc.  )  The  latter  are 
characteristic  of  a  system  reaching  a  coexistence  boundary  at 
other  conditions  (bubble  points  or  dew  points)  than  the 
critical.  Figure  5.11  displays  typical  experimental  results 
for  a  near-critical-density  gas  mixture.  Curve  1  represents  a 
run  for  which  the  bath  temperature  was  ramped  downward  from 
-30  to  -12  °C.  Curve  2  shows  data  for  an  upward  ramp  from  -10 
to  -22  °C.  In  both  cases  the  tracked  peak  (first  radial  mode 
resonance  f reguency)  was  lost  over  a  small  temperature 
interval  in  the  critical  vicinity  due  to  the  extremely  fast 
change  of  sonic  speed  and  drop  in  amplitude.  Curve  3  is  like 
curve  2  except  that  the  temperature  of  the  well-insulated 
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system  was  allowed  to  raise  naturally  from  15  °C  to  room 
temperature.  This  procedure  was  similar  to  the  one  utilized  in 
the  experiments  with  carbon  dioxide  (and  also  with  ethane  as 
will  be  seen  later)  except  that  the  run  went  from  low  to  high 
temperature  as  the  critical  temperature  of  this  particular 
mixture  is  lower  than  room  temperature.  Again,  the  peak  is 
easier  to  track  if  the  system  temperature  changes  very  slowly. 
The  speed  of  sound  is  very  sensitive  to  temperature  change  in 
the  critical  region.  Note  that  a  sharp  turning  point  gives  a 
guite  accurate  value  of  the  transition  temperature  but  not 
necessarily  of  speed  of  sound.  Figure  5.12  is  similar  to 
figure  5.11  but  illustrates  how  other  modes  of  vibration 
reveal  the  same  critical  point  temperature.  In  principle  every 
normal  mode  resonance  exhibits  a  discontinuity  at  the  critical 
point,  so  there  are  very  many  possible  indicators  which  are 
egually  capable  of  locating  the  critical  point. 

Figure  5.13  shows  the  behavior  of  pressure  versus 
temperature  for  several  isochoric  runs  representing  different 
fluid  densities.  The  upward  direction  of  the  breaks  in  these 
curves  shows  them  each  to  result  from  approaching  the  phase 
boundary  on  cooling  from  a  supercritical  density.  The  breaks 
mark  individual  bubble  points  on  the  phase  envelope. 
Approaches  from  subcritical  densities  would  break  downward  at 
the  corresponding  dewpoints,  and  the  critical  isochore  would 
exhibit  no  break  at  all.  These  facts  are  helpful  in  searching 
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for  the  critical  point  of  a  system  when  it  is  not  previously 
known.  For  example  the  starting  pressure  at  30  °C  for  the 
critical  isochore  could  be  approximately  calculated  from  the 
information  in  figure  5.13  in  the  following  fashion.  For  each 
isochore  the  two  branches,  supercritical  vapor  above  the 
bubble  point  and  two  phase  equilibrium  below  it,  were  fit  to 
straight  lines  by  linear  regression.  The  angle  0  between  the 
two  branches  is  directly  related  to  the  slopes  of  these  lines. 
A  plot  of  initial  charging  pressure  versus  0  indicates  by 
extrapolation  the  approximate  critical  density  charging 
pressure  for  which  0=0.  The  results  for  the  data  of  figure 
5.13  are  tabulated  in  table  5.3,  and  the  plot  of  charging 
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Figure  5.13.  A  sequence  of  supercritical  pressure  versus 
temperature  isochores  of  the  carbon  dioxide-ethane 
mixture. 
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pressure  versus  6  is  shown  in  figure  5.14.  This  technique 
revealed  the  magnitude  of  the  charging  pressure  to  be  used 
when  data  were  to  be  collected  in  search  of  the  critical 
point. 

Table  5.3.  Angles  between  pressure-temperature  lines  before 
and  after  phase  boundary  of  supercritical  density  fluid  of 
carbon  dioxide  and  ethane. 


File 

# 

Slope  of  line 

before  phase 

boundary 

Slope  of  line 

after  phase 

boundary 

Angle  between 

two  lines 

(degree) 

729 

34.637 

18.307 

14.733 

731 

34.074 

18.293 

14.479 

801 

33.562 

18.371 

14.091 

802 

32.842 

18.444 

13.594 

807 

27.882 

18.690 

10.356 

811 

25.149 

18.607 

7.992 

815 

24.701 

18.846 

7.191 

This  prediction  was  important  here  since  the  system  volume  is 
fixed.  The  only  way  to  change  the  system  density  was  to 
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Figure  5.14.  Relationship  between  the  angles  6  and  the 
charging  pressures  of  the  carbon  dioxide-ethane  mixture. 

release  some  of  the  material,  and  this  process  was  potentially 
troublesome  because  of  the  possibility  of  inducing  unwanted 
composition  changes.  Care  was  taken  to  release  gas  only  from 
the  system  in  the  single  phase  condition  and  then  only  after 
it  had  been  blended  thoroughly  by  continued  mixing.  Typically, 
when  the  system  was  close  to  the  critical  density,  only  small 
amounts  of  gas  were  released  for  each  successive  run.  The  drop 
in  charging  pressure  at  30  °C  was  typically  one  or  two  psia. 
For  this  work  the  total  number  of  runs  was  fifty. 

Table  5.4  lists  the  experimental  results  needed  to 
construct  phase  diagrams  of  this  system.  Figures  5.15  to  5.19 
show  various  types  of  phase  diagram  created  from  different 
combinations  of  the  parameters  in  table  5.4.  Figure  5.15  shows 
the  relationship  between  starting  pressure  at  30  °C  and  the 
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84 
temperature  of  phase  change.  Experimental  results  are 
represented  by  points.  The  smooth  curve  is  a  result  of  a  third 
order  polynomial  curve  fitting  procedure.  The  corresponding 
equation  is  as  follows: 

T  =  a  +  b*p  +  c*p2  +  d*p3 
where  T  =  temperature,  Celsius 

p  =  pressure,  psia 

a  =  -8.137773  X  102 

b  =  2.020012 

c  =  -1.616095  X  lO'3 

d  =  4.260451  X  10"7 

r  =  correlation  coefficient 
=  0.99749. 
Figure  5.16  shows  the  relationship  between  starting  pressure 
and  the  speed  of  sound  at  the  phase  change.  A  sharp  change  in 
the  speed  of  sound  at  the  critical  point  is  evident  compared 
to  the  slow  change  of  temperature  shown  figure  5.15.  Figure 
5.17  shows  a  coexistence  curve.  Note  that  the  curve  is  so  slim 
that  the  bubble-point  line  almost  overlaps  the  dew-point  line. 
Recall  that  a  similar  plot  for  a  pure  fluid  yields  a  nearly 
straight  line  (see  figure  5.23  of  ethane.)  This  clearly 
reveals  a  property  of  this  mixture  as  being  nearly  azeotropic. 
It  behaves  much  as  if  it  were  a  pure  fluid.  The  literature 
value  of  the  azeotrope  composition  of  this  binary  mixture  is 
reported  to  be  0.255  mole  fraction  of  ethane(108).  This  is  very 
close  to  the  composition  of  the  mixture  studied  in  this  work 
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Figure  5.15.  Relationship  between  the  starting  pressures  and 
the  temperatures  of  phase  change  of  the  carbon  dioxide-ethane 
mixture. 
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Figure  5.16.  Relationship  between  the  starting  pressures  and 
the  speed  of  sound  at  a  phase  change  of  the  carbon  dioxide- 
ethane  mixture. 
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Figure  5.18.  The  sonic  speed  versus  temperature  of  the  carbon 
dioxide-ethane  mixture  for  an  isochore  near  critical  density. 


940 


920 


900 


o   880 

1 

£    860 


840 


820 


15     16     17     18     19     20     21 

Temperature.  Celsius 


Figure  5.19.  Pressure  and  temperature  behavior  of  the  carbon 
dioxide-ethane  mixture  charged  near  its  critical  density. 
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Figure  5.20.  Pressure  and  temperature  relationship  of  a  set  of 
several  isochores  of  the  carbon  dioxide-ethane  mixture. 
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90 
(0.256)  confirming  the  explanation.  Figures  5.18  and  5.19  show 
the  relationships  near  the  critical  density  between 
temperature  of  phase  change  and  speed  of  sound  and  the 
pressure  respectively.  As  expected  the  sonic  speed  (figure 
5.18)  shows  a  pronounced  change  in  the  critical  region  while 
the  pressure  (figure  5.19)  shows  a  smooth  passage  from  the 
one-phase  to  the  two-phase  region.  Figure  5.20  shows  a 
composite  plot  of  pressure  and  temperature  data  representing 
23  separate  isochoric  runs.  Those  loadings  with  supercritical 
densities  show  phase  change  by  upward  breaks  at  the 
corresponding  bubble  points.  Those  loadings  at  subcritical 
densities  show  downward  breaks  at  the  corresponding  dew 
points.  Loading  at  the  critical  density  leads  to  no  break  at 
all.  Table  5.5  compares  dew  point  pressures  observed  in  this 
work  with  those  available  in  the  literature.  Clearly,  the 
agreement  is  guite  good.  The  critical  temperature  for  this 
mixture  obtained  in  this  work  was  293.13410.007  K.  The  only 
available  literature  experimental  value  for  this  composition 
was  found  to  be  293.072  K  <96'.  These  values  differ  by  only 
0.021  percent. 


Ethane 


The  graphs  shown  in  figures  5.21  to  5.24  represent  the 
isochoric  behavior  of  the  system  charged  with  pure  ethane  at 
near  critical  density.  Figure  5.21  shows  a  typical 
relationship  between  a  resonance  freguency  and  the  system 
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Figure   5.21.     Resonance    frequency     (Is    mode)     and    temperature 
relationship   of   ethane   near   the  critical   density. 
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Figure  5.22.   Resonance  frequency   (Is  mode)   and  pressure 
relationship  of  ethane  near  the  critical  density. 
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Figure  5.23.  Pressure  and  temperature  relationship  of  Is  mode 
resonance  frequency  of  ethane  near  the  critical  density. 


230 


Q 

o 

00 


o 

Q 


a. 

oo 


220- 

210- 

200 

190 

180 

170  H 

160 

150 

140 


130 


28  29  30  31  32  33  34  35  36  37  38 
TEMPERATURE,  CELSIUS 


Figure  5.24.   Speed  of  sound  and  temperature  of  Is  mode 
resonance  frequency  of  ethane  near  the  critical  density. 


93 
temperature.  Figure  5.22  is  similar  to  figure  5.21  but  with 
pressure  plotted  instead  of  temperature.  These  two  figures 
have  the  same  characteristics  implying  a  linear  relationship 
between  temperature  and  pressure  as  expected  for  this 
condition.  This  is  even  more  evident  in  figure  5.23.   Figure 


Table  5.6.  Sixth  order  polynomial  fit  statistics  of  GRAPHER 
program  applied  to  experimental  data  points  in  figure  5.25. 


Degree 

Polynomial 
Coefficients 

Sums  of  Squares 
of  Residuals 

0 

2.15029E+007 

0.952964 

1 

-166034 

0.909732 

2 

533.701 

0.0191329 

3 

-0.914178 

0.00389193 

4 

0.000880105 

0.0019246 

5 

-4.51543E-007 

0.00192265 

6 

9.64558E-011 

0.00192195 

Total  points  =23 

Points  in  fit  interval  =  23 

5.24  shows  plot  of  the  speed  of  sound  as  a  function  of 
temperature.  It  is  similar  to  figure  5.21  since  the  speed  of 
sound  is  directly  proportional  to  the  resonance  frequency. 
Figure  5.25  shows  the  relationship  between  the  initial 
charging  pressure  and  the  coexistence  (liquid-vapor 
equilibrium)  temperature.  Experimental  data  points  were  fitted 
to  a  polynomial  of  the  sixth  order.  Its  coefficients,  along 
with  some  statistical  values  generated  from  the  data  by  GRAPHER 
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Figure  5.25.  Phase  diagram  of  pure  ethane  showing  temperature 
of  phase  changes  as  a  function  of  starting  pressure. 


software,  are  listed  in  table  5.6.  Figure  5.26  shows  a 
collection  of  pressure  and  temperature  data  of  23  isochores 
which  cover  a  broad  range  of  densities  from  supercritical  to 
subcritical  values.  These  curves  coincide  within  experimental 
uncertainty  in  the  liquid-vapor  equilibrium  region.  Figure 
5.27  shows  a  trajectory  plot  of  a  collection  of  23  pressure- 
temperature-speed  of  sound  data  sets  of  different  densities. 
It  clearly  shows  the  pronounced  trend  of  decreasing  sonic 
speed  as  the  system  approaches  critical  conditions.  Table  5.7 
is  a  compilation  of  all  necessary  experimental  data  to 
construct  the  phase  diagram  in  different  views  as  shown  in  the 
above  listed  figures.  Tables  5.8  gives  the  comparison  between 
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Figure  5.26.  Phase  diagram  of  ethane  showing  pressure  and 
temperature  behavior  of  23  different  isochores. 
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Figure  5.27.  Trajectory  plot  of  temperature,  pressure,  and 
sonic  speed  of  ethane  at  different  densities.  The  one  nearest 
the  critical  density  reaches  the  lowest  sonic  speed. 


97 
the  present  experimental  bubble  point  pressures  and  available 
literature  values.  These  include  the  experimental  values 
measured  in  the  compressibility  apparatus  of  Douslin  and 
Harrison'109'  and  the  values  predicted  by  DDMIX(5S),  a 
correlation  model  developed  at  the  National  Institute  of 
Standards  and  Technology  (NIST,  formerly  NBS).  Table  5.9 
compares  the  present  vapor  pressures  with  those  of  the 
predictive  equation  proposed  by  Sychev  et.  al.(110)  This 
equation  derived  from  a  compilation  of  vapor  pressures  of 
ethane  by  Goodwin  et.  al.(111)  The  Sychev  equation  is: 

In  ps  =  a  +  bx  +  cu  +  du2  +  eu3  +  fu(l-u)€ 
where  x(T)  =  [  l-(Ttr/T)  ]  /  [  1-  (Ttr/Tcr)  ] 
u(T)  =  [T-Ttr]/[Tcr-TtJ 
Ttr  =  Temperature  at  the  triple  point 

=  90.348  ±  0.010  K 
Tcr  =  Critical  Temperature 
=  305.33  ±  0.02  K 
a  =  -11.38996 
b   =  18.84523 
c   =  -7.635425 
d  =  5.428443 
e  =  -1.362327 
f   =  0.7692493 
e   =  1.30. 
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Table  5.9.  Comparison  of  experimental  vapor  pressures  and 
those  calculated  with  an  equation  of  state  of  ethane  proposed 
by  Sychev.(110) 


Experimental 

Predicted 

pressure 

psia 

Percent 
difference 

T  (K) 

P  (psia) 

304.72 

697.65 

697.36 

0.0410 

304.95 

700.83 

700.79 

0.0053 

305.14 

703.71 

703.66 

0.0062 

305.16 

704.30 

703.96 

0.0481 

305.19 

704.70 

704.45 

0.0360 

305.24 

705.40 

705.16 

0.0338 

305.27 

706.09 

705.55 

0.0764 

305.30 

706.39 

706.05 

0.0475 

305.30 

706.59 

706.08 

0.0713 

305.31 

706.79 

706.27 

0.0735 

305.32 

705.89 

706.42 

-0.0747 

305.32 

707.08 

706.45 

0.0896 

305.33 

707.08 

706.47 

0.0874 

305.35 

707.78 

N/A* 

N/A 

305.36 

707.78 

N/A* 

N/A 

305.31 

706.49 

706.19 

0.0421 

305.24 

705.40 

705.23 

0.0231 

305.18 

704.50 

704.30 

0.0292 

305.10 

703.11 

702.99 

0.0176 

304.96 

700.53 

700.98 

-0.0648 

304.83 

698.54 

699.05 

-0.0726 

304.62 

696.16 

695.99 

0.0242 

*  Equation  diverges  at  this  temperature. 
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The  predictive  and  experimental  pressures  are  in  good 
agreement.  Notice  that  this  eguation  is  indeterminate  when 
temperature  is  greater  than  305.33  K  since  the  u(T)  term  is 
more  than  unity.  Table  5.10  chronologically  collects  the 
literature  values  of  the  critical  parameters  of  ethane  and  the 
value  obtained  in  this  work. 

Analysis  of  Ethane  Data 
As  mentioned  earlier  in  the  section  on  experimental 
procedure,  the  ethane  experiment  was  performed  differently 
from  those  on  carbon  dioxide  and  the  carbon  dioxide-ethane 
mixture.  This  resulted  from  advances  in  the  experimental 
technique  as  the  research  progressed.  In  the  case  of  ethane 
all  data  points  along  the  entire  frequency  sweep  range  were 
saved  on  memory  discs.  These  data  can,  in  turn,  be  examined 
and  analyzed  later  to  obtain  more  understanding  of  the 
behavior  of  the  gas  under  acoustic  excitation.  The  analysis, 
however,  is  still  very  difficult  and  time  consuming  due  in 
part  to  the  large  number  of  data  points.  Graphic  methods  are 
appropriate  for  this  purpose.  For  example,  in  the  ethane 
experiment  for  each  isochoric  density  there  are  approximately 
400  data  files  each  containing  120  sets  of  four  parameters: 
temperature,  pressure,  resonance  frequency,  and  amplitude.  For 
two  dimensional  observations  one  would  have  to  pick  two  out  of 
four  parameters  to  plot  each  data  set  for  the  total  of  120X400 
plots.  The  practicability  of  doing  this  is  limited  by  many 
factors  such  as  the  available  graphics  capability,  computer 
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Figure  5.28.  Resonance  frequency  and  temperature  relationship 
revealing  some  phenomena  observed  in  the  experiment  on  pure 
ethane. 

memory  size  and  resolution  of  the  display  monitor  etc.  At 

present  only  limited  observations  were  practical. 

Figure  5.28  is  used  in  conjunction  with  the  following 
observations.  Note  that  each  curve  in  this  figure  does  not 
necessarily  represent  a  system  having  the  same  density  as  the 
others.  They  are  considered  together  for  the  purpose  of  the 
following  comparison. 
Observation  1 

At  temperatures  above  the  critical  temperature  the  first 
radial  mode  of  vibration  (the  Is  peak)  and  the  first  nonradial 
mode  of  vibration  (the  If  peak)  were  in  close  proximity  (about 
15-20  hz  apart).  Consequently,  the  If  peak  tended  to  interfere 
with  tracking  of  the  Is  peak.  The  data  acquisition  program 
selected  the  signal  of  highest  amplitude  in  the  observation 
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window,  and  this  was  prone  to  shift  between  these  two  on 
occasion.  This  situation  caused  an  unsmooth  experimental  curve 
as  shown  in  the  curve  labeled  observation  1  of  figure  5.28. 
Observation  2 

On  cooling  below  the  critical  temperature  in  some  runs 
the  tracking  program  locked  onto  some  peak  other  than  the 
desired  Is  peak  and  showed  aberrant  behavior  on  continued 
cooling.  The  curve  labeled  observation  2  of  figure  5.28  shows 
a  behavior  of  this  sort.  No  further  investigation  was  done  to 
identify  the  spurious  signal. 
Observation  3 

If  the  situation  in  observation  2  did  not  occur,  the  peak 
emerging  after  the  critical  condition  is  the  Is  peak3.  At 
still  lower  temperatures,  however  this  is  lost  intermitentally 
to  a  companion  peak  which  grows  in  amplitude  beyond  that  of 
the  Is  peak.  The  partial  history  of  the  scanned  signal  is 
depicted  in  the  plots  of  figures  5.29  and  5.30  which  show  the 
splitting  of  the  tracked  peak.  The  scans  are  shown  in  proper 
time  order  but  those  for  which  the  peak  is  single  are  mainly 
omitted.  Curves  4  and  5  of  figure  5.29  show  the  onset  of 
splitting.  As  the  system  moved  away  from  the  critical  point 


3  Based  on  the  technigue  used  in  this  work  one  cannot  say 
confidently  what  the  emerged  peaks  are.  The  only  sure  way  to 
solve  this  is  to  collect  all  peaks  in  the  entire  freguency 
range  (from  0  to  10,000  khz)  and  observe  the  relationships 
among  them.  This,  however,  was  beyond  the  scope  of  this  work. 
Other  work  in  this  laboratory,  however,  is  directed  to  this 
problem. 
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Figure  5.29.  Movement  of  tracked  peak  and  its  neighbor  below 
the  critical  temperature.  These  curves  are  numbered  in  the 
proper  time  sequence  but  many  curves  showing  only  one  strong 
peak  have  been  omitted. 
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Figure  5.29.  Continued  movement  of  tracked  peak  and  its 
neighbor  below  the  critical  temperature. 


110 
the  separation  of  these  two  peaks  increased  (curves  6  through 
11  of  figure  5.29.)  The  right  hand  peak  (peak  one)  was  bigger 
and  moved  faster  than  the  left  peak  (peak  two.)  It  moved  so 
fast  that  peak  two  was  eventually  out  of  the  sweep  range  (see 
curve  12  of  figure  5.29.)  This  single  peak  remained  strong 
until  the  temperature  was  around  31  °C.  At  that  temperature 
the  other  peak  (possibly  peak  two)  reappeared  in  the  sweep 
range  (curves  15  and  16  of  figure  5.29.)  It  was  at  this  time 
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Figure  5.30.  Possible  routes  of  movement  of  two  neighbor  peaks 
(peak  one  and  peak  two)  below  the  critical  temperature. 
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stronger  than  peak  one  leading  to  switching  of  the  tracked 
peak  from  peak  one  to  peak  two.  The  window  then  shifted  from 
peak  one  which  became  lost.  The  peak  in  curve  16  is  peak  one 
but  that  of  curve  17  is  peak  two.  This  is  apparently  so  since 
the  peak  in  curve  17  is  lower  in  frequency  than  the  peak  in 
curve  16.)  The  reverse  of  this  phenomenon  happened  around  29.5 
°C  where  peak  one  re-entered  the  window  as  the  dominant 
signal.  Finally,  these  two  peaks  merged  together  at  around  29 
°C.  Data  were  not  collected  below  about  29  °C,  so  the 
subsequent  behavior  of  these  peaks  was  not  observed.  This  type 
of  behavior  is  interesting  and  unexpected,  but  it  is  not 
central  to  the  objective  of  the  present  research.  It  is 
probably  associated  with  the  fact  that  below  the  critical 
temperature  there  are  two  phases  present  and  their  relative 
proportions  in  the  resonator  cavity  change  according  to  the 
effects  of  the  magnetically  driven  circulating  pump. 

The  plausible  routes  of  these  two  peaks  based  on  the 
above  situations  are  depicted  in  figure  5.30.  A  full 
interpretation  of  these  phenomena  would  require  further 
experimentation  to  characterize  the  contents  of  the  cavity  in 
the  two  phase  region.  In  the  early  stage  of  phase  separation 
the  bubbles  or  dew  drops  are  believed  to  be  small  and 
entrained  in  the  dominant  phase  as  foam  or  fog,  so  the 
location  of  the  phase  transition  is  not  obscured  by  massive 
phase  separations  which  occur  well  away  from  the  phase 
boundary. 
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Observation  4 

An  attempt  at  running  the  experiment  with  ascending 
temperature  was  made.  The  same  uncertainties  about  phase 
distribution  obscures  the  identify  of  the  peak  selected  for 
tracking.  The  possibility  of  selecting  a  peak  which  evolves 
irregularly  or  fades  into  the  background  exists.  The  curve 
labeled  observation  4  in  figure  5.28  shows  this  type  of 
observation. 

Curve  Fitting 

Over  a  small  temperature  interval  in  the  critical 
vicinity  the  task  of  single  peak  tracking  is  guite  severe,  and 
the  peak  is  often  lost.  This  is  due  to  two  main  effects. 
First,  in  terms  of  dynamics,  the  rate  of  decrease  in  resonance 
freguency,  hence  speed  of  sound,  with  respect  to  temperature 
is  extremely  high  approaching  infinity,  i.e.  the  speed  of 
sound  exhibits  a  mathematical  singularity  at  the  critical 
point.  Second,  the  amplitude  of  the  tracked  peak  decreases 
dramatically  toward  zero  or  at  least  below  the  noise  level  of 
the  system.  As  a  conseguence,  an  ambiguity  of  the  critical 
point  arises.  In  order  to  solve  this  problem  a  mathematical 
model  describing  the  relationship  between  the  speed  of  sound 
and  temperature  near  the  critical  was  applied.  This  model  was 
originally  derived  by  Stell  et.  al.(133)  to  fit  data  of  latent 
heat  of  vaporization  and  temperature  and  later  modified  to 
satisfy  a  corresponding  states  principle  by  Sivaraman  et. 
al#  (134,135)  and  appiied  to  speed  of  sound  data  by  the  same 
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Fxgure  5.31.  Experimental  curve  of  ethane  showing  dynamic 
behavior  of  speed  of  sound  near  critical  point. 

author.  It  has  a  root  from  rigorous  mathematical  derivations 
by  Wegner(136)  and  further  by  Green  and  Ley-koo(137'138)  based  on 
renormalization  group  theory. 
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(R)(TC) 

M 


A^  +  A2T;+A  ♦  A3Tl-a+p  ♦  A4Tr  +  A5T*   ♦  A6T^ 


(5) 


where   C   =  Speed  of  sound,  m/sec 
R  =  Gas  constant,  J/sec 
Tc  =  Critical  Temperature,  K 
M  =  Molecular  mass,  Kg/mole 
A±  =  Eguation  coefficients 
Tr  =  Reduced  Temperature 
a,    (S,    A  =  Critical  exponents. 
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To  fit  this  non  linear  equation  to  experimental  data  of 
speed  of  sound  and  temperature,  eleven  system-dependent 
coefficients  must  first  be  searched.  This  task  was 
accomplished  by  the  simplex  mathematical  optimization 
technique.  (A  computer  program  written  by  the  author  for  this 
purpose  is  given  in  appendix  C.)  In  general  the  experimental 
data  of  speed  of  sound  and  corresponding  system  temperature 
were  divided  into  two  parts  according  to  temperature.  The 
first  one  is  all  data  at  temperatures  above  the  critical 
temperature  (roughly  estimated  from  the  literature  value.)  The 
second  one  is  all  data  at  temperature  below  the  critical 
temperature.  Each  data  set  was  independently  fitted  to  the 
above  equation  resulting  in  a  set  of  eleven  coefficients. 


i  x  n 

ZjU 

8  220  - 

^S. 

<21  0  - 

^v 

E  200  - 

v 

o  1  9  0  - 

\         *+*****" 

8  180  - 

\    ^^ 

£  170  - 
o  1  60  - 

\/ 

S  1  50  - 
1 1 1 

f 

&   140  - 
1  z  n 

* 

B 

2 

3293031   3233343536373 

TEMPERATURE,     CELSIUS 

FITTING  CURVE                 +   EXPERIMENTAL     POINTS 

(PLOTTED  EVERY      4      POINTS) 

and  curve  generated  by  renormalization  group  theory  equation, 
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Temperature  and  speed  of  sound  at  the  intersection  of  these 
two  smooth  curves  were  treated  as  turning-point  parameters  of 
the  corresponding  experimental  data  set.  Table  5.11  shows 
results  of  searching  coefficients  of  experimental  data  shown 
in  figure  5.31  by  the  simplex  optimization  program.  Figure 
5.32  compares  the  experimental  data  and  fitting  curve 
generated  by  this  method  for  ethane  at  the  critical  density. 
Note  that  in  this  figure  only  one-forth  of  the  experimental 
points  were  plotted  to  preserve  picture  clarity.  Clearly,  this 
eguation  represents  very  well  the  behavior  of  the  speed  of 
sound  in  the  vicinity  of  the  critical  point.  Furthermore,  it 
also  yields  an  egually  satisfactory  result  for  the  mixture  of 
carbon  dioxide  and  ethane  described  in  the  previous  section. 


Table  5.11.  Coefficients  Obtained 
Method. 


by  Simplex  Optimization 


Coefficients 

High  Temperature 
Side 

Low  Temperature 
Side 

* 

11.95399178 

4.4228669 

A, 

-897.98691534 

360.1419335 

*■ 

-3400.8303361 

2448.9833655 

A< 

2976.6473117 

-1666.3105786 

A< 

7325.7816766 

-7444.4800181 

A« 

-56033.5889588 

52720.3509156 

T 

305.3979294 

305.9057717 

a 

0.1250001 

0.1250000 

P 

0.3333333 

0.3333333 

A 

0.5000000 

0.5000000 
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Three  Dimensional  Phase  Diagram 
Conventionally  the  phase  diagram  in  three  dimensions  is 
represented  by  using  parameters  of  temperature,  pressure,  and 
molar  volume  leading  to  three  well-known  two  dimensional 
projections:  isotherm,  isobar,  and  isochore  curves.  Figure  3.1 
illustrates  an  example  of  this  type  of  plot  for  a  typical  pure 
fluid.  The  critical  point  is  located  at  the  top  of  a 
coexistence  curve.  It  can  be  found  by  a  null  change  of  slope 
in  the  critical  isotherm  or  isobar  or  by  a  smooth  curve 
without  break  in  the  critical  isochore.  These  are  among  the 
traditional  indicators,  but  their  response  to  the  critical 
point  is  not  highly  sensitive. 

For  more  accurate  measurements  a  more  sensitive  indicator 
is  needed.  The  sonic  speed  or  cavity  resonance  frequency 
serves  nicely  for  this  purpose  due  to  their  dramatic  change  at 
the  critical  point.  Replacing  molar  volume  by  speed  of  sound 
in  the  conventional  phase  diagram  plot  leads  to  a  pyramid 
shape  curve  as  shown  in  figure  5.33  for  a  mixture  of  carbon 
dioxide  and  ethane.  Exploring  in  a  much  narrower  range  of 
temperature  around  the  critical  point  leads  to  an  outstanding 
tornado  shape  curve  as  shown  in  figure  5.34  for  ethane.  These 
two  curves  were  constructed  from  experimental  data.  The 
critical  point  is  apparently  easily  located  as  an  extrumum 
point  on  these  surfaces.  It  corresponds  to  the  point  of 
minimum  speed  of  sound,  which  theoretically  speaking,  should 
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be  zero.  In  that  case  P-T-C  phase  diagram  curve  would  have 
much  more  prominent  tornado-in-action  shape. 
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Figure  5.34.  Three  dimensional  phase  diagram  of  ethane, 


CHAPTER  6 
CONCLUSION 


A  new  approach  using  an  acoustic  resonance  technique  for 
phase  boundary  detection  may  be  summarized  as  follows: 

1.  Measure  several  resonance  frequencies  of  standing 
acoustic  waves  excited  inside  a  fluid  filled  spherical 
acoustic  cavity. 

2.  Identify  each  resonance  frequency. 

3.  Track  an  assigned  resonance  frequency  as  system 
parameters  such  as  temperature,  pressure,  or  density  are 
varied. 

4 .  Locate  phase  equilibria  by  abrupt  changes  in  resonance 
frequency  (hence  speed  of  sound)  as  the  system  state  traverses 
a  phase  boundary. 

5.  Locate  critical  points  as  minima  on  (P,T,C)  or  (p,T,C) 
diagrams . 

The  principal  results  of  this  work  are  tabulated  and 
compared  with  the  available  literature  values  in  table  6.1. 
Apparently,  the  present  results  are  in  good  agreement  with 
established  literature  values,  proving  the  viability  of  this 
new  technique.  The  largest  discrepancy  is  in  the  system 
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consisting  of  a  C02-C2H6  mixture.  This  is  not  surprising 
inasmuch  as  the  literature  value  was  obtained  by  extrapolation 
of  the  reported  experimental  data  and  thus  is  subject  to 
greater  uncertainty.  This  comparison  is  intended  principally 
to  show  the  capability  of  the  acoustic  method  to  detect  the 
critical  point  of  the  more  complicated  system  and  should  be 
judged  more  on  the  precision  than  the  accuracy  of  the  outcome 
since  the  basis  for  judging  the  latter  is  not  sufficiently 
strong. 

Generally  speaking,  considering  the  results  for  pure 
gases,  even  though  the  values  of  this  work  are  not  all  within 
the  reported  uncertainty  of  the  literature  values,  the  work 
should  not  be  judged  a  failure.  In  the  first  place  the  maximum 
disparity  between  the  different  methods  is  only  1  mK,  and  an 
increase  in  the  estimated  limits  of  uncertainty  of  only  ±1  mK 
would  bring  the  divers  results  into  accord.  There  are  reasons, 
however,  to  suspect  that  the  present  results  may  be  more 
reliable.  Consider,  for  example,  the  conventional  method  of 
detecting  the  critical  point.   This  involves  the  direct 
observation  of  the  transition  from  two  phases  to  one,  with  the 
critical  state  taken  to  be  the  single  condition  for  which  the 
meniscus  disappears  when  the  volumes  of  liguid  and  vapor  are 
just  the  same.  Ideally,  critical  points  would  be  determined  by 
the  direct  observation  of  this  phenomena  at  the  exact  middle 
of  the  cylindrical   sample  cell.   Realistically,   such  a 
measurement   is   experimentally   impracticable*96'   and   the 
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appearance  of  the  meniscus  elsewhere  is  used  as  a  to  determine 
the  critical  parameters.  In  practice  these  are  determined  by 
interpolation  of  data  taken  when  the  meniscus  appears  just 
above  and  just  below  the  middle  of  the  cell.  This  method  is 
then  subject  to  an  heuristic  effect.  Besides,   near  the 
critical  condition  the  phenomena  of  critical  opalescence 
(light  scattering  effect)  may  obscure  the  observation  to  some 
extent.  In  the  acoustic  method,  the  instrument  is  generally 
capable  of  measuring  freguencies  with  an  accuracy  of  at  least 
1  mHz.  In  the  course  of  these  experiments  data  were  taken  over 
a  4  K  range  around  the  critical  temperature  over  a  period  of 
twelve  hours.  The  recorded  data  show  that  nearest  the  observed 
critical  state  the  temperature  gradient  of  resonance  f reguency 
exceeds  2000  Hz/K  on  both  sides  of  the  freguency  minimum  for 
which  the  corresponding  temperature  was  taken  as  a  critical 
temperature.  A  remarkably  abrupt  change  in  freguency  is 
evident.  The  sharpness  of  this  turning  point  together  with  the 
stated  precision  of  freguency  measurements  illustrates  how  the 
present  technigue  can  locate  the  extremum  to  a  resolution  of 
about  one  part  per  million  in  the  temperature.  Demonstrating 
this  fact  was  indeed  the  principal  goal  of  this  research.  That 
the  acoustic  resonance  method  is  vastly  more  sensitive  (ca. 
three  orders  of  magnitude)  than  the  visual  method  is  far  more 
important  than  whether  of  not  the  accuracy  of  temperature 
measurement  at  a  small  poorly  funded  university  laboratory  is 
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comparable  to  that  at  the  major  standards  laboratory  of  the 
nation.  Given  equivalency  of  thermometric  accuracy,  the  new 
method  of  critical  point  determination  is  clearly  superior  and 
is  a  good  candidate  to  become  the  method  of  choice  in  future 
studies.  Even  so  improvements  are  possible  and  some  of  these 
are  being  incorporated  in  the  next  generation  version  of  this 
apparatus.  The  major  change  will  be  the  inclusion  of  a  volume 
change  capability  so  that  data  can  be  taken  isothermally  or 
isobarically  as  well  as  isochorically. 

The  spherical  resonator  is  a  remarkably  accurate  and 
convenient   tool   for   the   measurement   of   thermophysical 
properties  of  fluids.  The  speed  of  sound  in  a  gas  of  interest 
can  be  measured  with  high  accuracy  merely  by  measuring  the 
frequencies  of  the  radial  normal  modes  of  vibration.  A 
relative  measurement  against  a  reference  gas  such  as  argon 
further  simplified  the  procedure.  This  technique  offers  a 
proven  alternative  experimental  approach  to  valuable  physical 
properties    including   phase   boundaries,    thermodynamic 
properties  and  equation  of  state  parameters,  all  of  which  are 
useful  to  both  pure  and  applied  science.  Applications  are 
numerous  and  only  now  beginning  to  be  explored.  One  important 
application  is  indirectly  determining  the  intermolecular 
potentials  of  fluids  from  acoustic  virial  coefficients  via  an 
inversion  technique  without  using  the  parameter  optimization 
techniques   for  model  potentials   such  as  Lennard- Jones . 
Accurate  sonic  speed  data  as  a  function  of  temperature  of  both 
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pure  fluids  and  mixtures  would  be  useful  source  of  this 
information.  The  properties  of  bulk  matter  derive  ultimately 
from  the  interactions  between  its  molecules,  and  information 
about  these  interactions  supports  the  development  of  useful 
correlations  and  of  statistical  mechanics  in  general. 

It  is  apparent  from  this  work  and  others  that  the  natural 
vibrations  of  fluid  masses  convey  a  vast  storehouse  of  encoded 
information.  Learning  to  acquire  and  interpret  that 
information  is  the  long  range  goal  of  this  laboratory.4 


Some  work  has  already  begun.  See  reference  139 


APPENDIX  A 
CHARGING  PRESSURE  CALCULATIONS 


Since  in  this  work  data  were  collected  in  the  isochoric 
mode,  the  density  of  the  system  was  constant  throughout  each 
run  as  well.  To  locate  the  critical  without  further  addition 
of  gas  required  that  the  initial  runs  be  made  on  samples  at 
supercritical  densities.  Density  was  not  directly  measured  but 
was  infered  from  the  AGA8  equation  of  state(5);  therefore  the 
starting  conditions  were  specified  by  composition, 
temperature,  and  pressure.  The  procedures  followed: 
Pure  Fluid  (Carbon  dioxide  or  Ethane) 

1.  Obtain  literature  values  of  the  critical  point 
parameters  (Tc,  Pc,  pc.) 

2 .  Use  the  AGA8  equation  to  calculate  a  series  of 
densities  corresponding  to  varying  pressures  in  the  vicinity 
of  the  critical  pressure  at  a  charging  temperature  which  must 
be  higher  than  the  critical  temperature  by  about  five  to  seven 
degrees  celsius.  This  set  of  densities  should  span  the 
critical  density. 

3.  Plot  the  calculated  densities  as  a  function  of  the 
corresponding  charging  pressures,  then  apply  a  suitable  curve 
fit  regression,  generally  polynomial  of  the  nth  order,  to 
these  data. 
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4.  Charging  pressure  corresponding  to  literature  critical 
density  is  now  calculated  from  the  best  curve  fitting  equation 
found  in  step  3.  This  value  plus  about  10  psia  is  taken  as  a 
charging  pressure  of  fluid  at  the  starting  temperature.  The 
extra  10  psia  is  added  as  a  safety  factor  due  to  the 
possibility  of  error  in  the  literature  critical  density. 
Fluid  Mixture 

1.  Decide  on  the  composition  of  carbon  dioxide-ethane 
mixture  to  be  studied. 

2 .  Approximate  charging  pressure  at  starting  temperature 
from  closest  literature  critical  point  parameters  available. 

3.  Calculate  the  density  of  the  gas  mixture  corresponding 
to  the  charging  pressure  estimated  in  step  2  at  the  starting 
temperature  and  mole  fractions  of  carbon  dioxide  and  ethane 
using  the  AGA8  equation. 

4.  Calculate  the  density  of  ethane  from  the  density  of 
the  gas  mixture  obtained  from  step  3. 

5.  Follow  procedures  in  step  2  to  4  for  a  pure  fluid.  The 
calculated  density  obtained  in  step  3  for  this  procedure  is 
taken  as  the  literature  critical  density  called  for  in  step  4 
of  pure  fluid.  The  resulting  pressure  is  taken  as  the  target 
charging  pressure  of  ethane. 

6.  Ethane  was  charged  first  to  the  pressure  calculated  in 
the  previous  step.  Carbon  dioxide  was  subsequently  introduced 
to  the  resonator  while  mixing  until  the  desired  total  pressure 
was  obtained  (estimated  value  in  step  2.) 


APPENDIX  B 
QUANTITATIVE  ANALYSIS  OF  GAS  MIXTURE 


Reported  values  of  critical  point  parameters  for  gas 
mixtures  are  meaningless  unless  accompanied  by  their 
quantitative  compositions.  For  this  particular  mixture  of 
carbon  dioxide  and  ethane  two  methods  were  used  for 
quantitative  analysis:  a  gravimetric  method  for  carbon  dioxide 
and  a  gas  chromatography  method  for  ethane.  Carbon  dioxide  was 
not  analyzed  by  the  latter  method  due  to  an  unavailability  of 
a  suitable  detector  during  the  course  of  analysis.  The 
assumption  used  through  out  this  analysis  was  that  there  were 
no  significant  impurities  since  the  mixture  was  carefully 
prepared  in  this  laboratory  from  highly  pure  gases. 
Gravimetric  Method 

Carbon  dioxide  reacts  with  base  to  form  carbonate  salt 
which  readily  dissolves  in  water  as  shown  in  the  following 
chemical  equation: 

C02  +  2  OH"  =  C032"  +  H20 
The  difference  between  weights  of  reaction  vessels  before  and 
after  carefully  passing  carbon  dioxide  is  the  weight  of  the 
gas.  The  weight  of  ethane  was  calculated  by  subtracting  the 
measured  weight  of  carbon  dioxide  from  the  weight  of  the 
sample . 
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Figure  B.l.  Block  diagram  of  the  gravimetric  method. 

A  block  diagram  of  the  experimental  setup  is  shown  in 
figure  B.l.  A  gas  sampling  container  was  made  from  a  stainless 
steel  tube  provided  with  high  pressure  needle  valves  A  and  B. 
The  experiment  started  with  pre-conditioning  the  system  by 
flowing  argon  gas  to  the  reaction  vessels  (hexagonal  absorber 
array  5,  see  figure  B.2)  through  valves  4  and  5  for  one-half 
hour.  After  that  valves  4  and  5  were  closed  followed  by 
opening  valve  B.  The  pressure  of  gas  in  the  gas  sampling 
container  was  shown  on  gauge  3  (typical  pressure  was  around 
1200  psia.)  Then  valve  5  was  opened  carefully  to  bleed  gas 


The  hexagon  absorber  arrays  of  vertical  tubes  each 
having  a  wash  board  shape  were  designed  to  prolong  contact 
time  of  gas  with  absorbant  to  assure  completeness  of  the 
reaction.  Each  section  of  the  first  five  tubes  was  filled  with 
1M  NaOH.  The  last  tube  was  filled  with  drierite  (moister 
absorbent.) 
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Figure  B.2.  The  hexagon  reaction  vessel  for  a  gravimetric 
analysis. 
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Figure  B.3.  The  high  pressure  gas  sampling  container. 
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slowly  to  the  hexagon  reaction  vessels.  Opening  of  valve  5  was 
judged  from  the  rate  of  bubbles  rising  up  in  the  absorber 
tubes.  If  valve  5  was  opened  too  wide,  carbon  dioxide  gas 
would  not  have  time  to  react  efficiently  with  sodium  hydroxide 
in  the  absorber.  The  completeness  of  the  reaction  was  checked 
by  passing  exiting  gas  from  the  second  absorber  array  through 
a  lime  water,  saturated  Ca(OH)2  solution.  Turbidity  of  the 
solution  due  to  forming  of  CaC03  indicates  incompleteness  of 
the  reaction.  Once  bubbles  stopped  rising  up  (about  2-3  hours) 
argon  was  delivered  to  the  system  through  valve  A  to  flush  out 
all  sample  gas  remaining  in  the  lines  for  one-half  hour.  After 
that,  screw  clamps  of  the  entrance  and  exit  tygon  tubes  to  the 
absorber   array   were   tightly   closed.   The   vessels   were 
disconnected  from  the  system  and  weighed  on  an  analytical 
balance.  During  the  entire  course  of  these  experiments  the 
lime  solution  remained  clear.  The  total  absorption  time  was 
about  3-4  hours.  Note  that  before  analyzing  the  sample  gas  a 
blank  test  was  performed  by  continuously  flowing  argon  through 
the  system  for  3-4  hours.  No  significant  change  in  weight 
(comparable  with  the  uncertainty  of  the  balance,  1  mg.)  of  the 
absorbers  before  and  after  passing  argon  was  observed.  Despite 
its  tedious  procedure  this  method  worked  reasonably  well, 
provided  that  sufficient  care  was  exercised. 
Gas  Chromatography  Method 

Description  of  standard  procedures  for  performing  gas 
chromatography  can  be  found  elsewhere.  (140'141»  A  calibration 
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curve  for  ethane  was  made  first  using  pure  ethane.  Both 
calibration  and  sample  analysis  were  performed  on  the  same  day 
basis  to  minimize  errors.  Details  of  the  equipment  and 
optimized  conditions,  which  were  determined  by  trial  and  error 
methods,  are  summarized  as  follows: 

Gas  chromatograph:  Varian  model  3700 


Syringe:  Pressure-Lock  gas  syringe  series  A  with  side 
port.  Dynatech  Precision  Sampling  Corporation. 

Column:  Fused  silica  megabore  column  GSQ,  0.53  mm  i.d., 
30  m  long,  open  tubular.  J&W  Scientific. 

Optimized  conditions: 

Temperature (°C)  -  Injector:  100 

-  Flame  Ionization  Detector:  100 

-  Column  :  80 
Flow  rate(ml/min)  -  Helium  :  6 

-  Hydrogen  :  30 

-  Air  :  300 

Count  rate  on  the  integrator  :  12,000  cpm. 
Sample  results  of  calibrations  are  summarized  in  table  B.l. 
Figure  B.3  shows  the  calibration  curve  for  ethane. 


Table  B.l.  Data  for  calibration  curve  of  ethane. 
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APPENDIX  C 
SIMPLEX 

There  are  three  programs  for  searching  for  minimum  points 

of  speed  of  sound  (points  along  coexistence  curve).  The  first 

is  the  simplex  optimization  program  used  to  search  for  eleven 

coefficients  in  an  equation  (5).  The  second  program  is  used  to 

calculate  speed  of  sound  with  coefficients  found  in  program 

one.  The  last  program  finds  the  common  point  of  two  curves 

generated  from  programs  one  and  two. 

Program  I.  Simplex  Optimization 

C234567 

IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  C(4000) ,T(4000) ,B(20,20) 

DIMENSION  SSRSD(500) 

WRITE (*,*)  'INPUT  MOLECULAR  WEIGHT  IN  KG/MOLE' 

READ  (*,*)  W 

OPEN  (10, FILE= ' DATA . RAW ' , STATUS= ' OLD ' ) 

L=l 

20    READ(10,*,END=30)  T(L),C(L) 
L=L+1 
GOTO  20 
30    CLOSE (10) 
L=L-1 

WRITE (*,*)  'L=',L 
WRITE (*,*)  'W=',W 

OPEN  (40, FILE= ' INITIAL . DAT ' , STATUS= ' OLD ■ ) 
READ(40,*)  ((B(I,J),I=1,10),J=1,11) 
CALL  INITLCOOR  (B,SSRSD,L,T,C,W) 
OPEN  ( 5 , FILE= ' START . DAT ' , STATUS= ' NEW ' ) 
DO  80  J-1,11 
DO  90  1=1,10 
WRITE  (5,*)  B(I,J) 
90   CONTINUE 

WRITE  (5,*)  •SSRSD,,J,'=',SSRSD(J) 
WRITE  (5,*)  ' 
80   CONTINUE 
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CLOSE ( 5 ) 
50   IF  (SSRSD(l) .LT.1.0E-17)  THEN 
C  WRITE (*,*)  'FIN  CON  1' 

GOTO  100 
END  IF 

CALL  SHELL (B, SSRSD, WORST) 
CHECKER=CHECKER+SSRSD ( 1 ) 
KOUNT=KOUNT+l 
C       WRITE (*  *)  ' *************************KOUNT= '  KOUNT 
KKOUNT=KKOUNT+ 1 
IF  (KKOUNT.EQ.300)  THEN 
C       IF  (KKOUNT.EQ.200)  THEN 
KKOUNT=0 

CHECKER= ( CHECKER/300 . 0 ) -SSRSD ( 1 ) 
C         CHECKER= ( CHECKER/200 . 0 ) -SSRSD ( 1 ) 
IF  (CHECKER. LT.1.00000E-17)  THEN 
C         IF  (CHECKER.LT. 1.0E-16)  THEN 
GOTO  100 
END  IF 
CHECKER=0.0 
END  IF 

IF  ( SSRSD ( 11 ) .EQ. WORST)  THEN 
C        WRITE (*,*)  'FLAG  RAISED  FLAG  RAISED  FLAG  RAISED' 
Z=SSRSD(10) 
SSRSD(10)=SSRSD(11) 
SSRSD(11)=Z 
DO  140  J=l,10 
ZZ=B(J,10) 
B(J, 10)=B(J, 11) 
B(J,11)=ZZ 
140    CONTINUE 
END  IF 
C       WRITE (*,*)  'START  REFLEX' 

CALL  REFLEX  (SSRSD, B,L,T,C,W) 
C       WRITE (*,*)  'PASS  REFLEX' 

IF  (SSRSD(13) .LT.SSRSD(l) )  THEN 
CALL  EXPAND(B, SSRSD, C,T,L,W) 
C        WRITE (*,*)  'PASS  EXPAND' 
GOTO  50 
ELSE  IF  (SSRSD(13) .GT.SSRSD(l) .AND. 
CSSRSD(13) .LT.SSRSD(IO) )  THEN 
SSRSD(11)=SSRSD(13) 
DO  70  1=1,10 
B(I,U)«B(I,13) 
70     CONTINUE 
C        WRITE (*,*)  'PASS  BETWEEN' 
GOTO  50 
ELSE  IF  (SSRSD (13) .GT. SSRSD (10) )  THEN 
IF  (SSRSD(13) .LT.SSRSD(ll) )  THEN 
CALL  CNTR(B, SSRSD, C,T,L,W) 
C         WRITE (*,*)  'PASS  CNTR' 

ELSE  IF  (SSRSD(13) .GT.SSRSD(ll) )  THEN 
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C 
C 

c 


CALL  CNTW(B, SSRSD, C,T,L,W) 
WRITE (*,*)  'PASS  CNTW' 
END  IF 
GOTO  50 
END  IF 
CLOSE (40) 
100  OPEN  (60, FILE= ' COOR . OUT ' , STATUS= ' NEW ' ) 
DO  800  J=l,ll 
DO  900  1=1,10 
WRITE  (60,*)  B(I,J) 
900  CONTINUE 

WRITE  (60,*)  'SSRSD' ,J, '=' ,SSRSD(J) 
WRITE  (60,*)  ' 
800  CONTINUE 

WRITE (60,*)  'NUMBER  OF  ITERATIONS  =  ',KOUNT 

CLOSE (60) 

STOP 

END 


SUBROUTINE  INITLCOOR  ( B, SSRSD, L,T,C,W) 
IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  B(20,20),T(4000),C(4000),SSRSD(500),A(10) 
DO  10  1=1,11 

DO  20  J-1,10 

A(J)=B(J,I) 

20    CONTINUE 

A1=A(1) 

A2=A(2) 

A3=A(3) 

A4=A(4) 

A5=A(5) 

A6=A(6) 

TC=A(7) 

ALPHA=A(8) 

BETA=A ( 9 ) 

DELTA=A(10) 

CALL  RESIDUAL (A1,A2, A3, A4,A5,A6,TC, ALPHA, BETA 
C, DELTA, C,T,L, ST, W) 
SSRSD(I)=ST 
C        WRITE(*,*)  ,SSRSD',I,'=',SSRSD(I) 


C 

c 

c 


10    CONTINUE 
RETURN 

END 


SUBROUTINE  RESIDUAL ( Al ,A2 , A3, A4 ,A5 ,A6 ,TC, ALPHA, BETA 
C, DELTA, C,T,L, ST, W) 
IMPLICIT  REAL*8  (A-H,0-Z) 
REAL* 8  DABS,DSQRT 
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DIMENSION  C(4000) ,T(4000) 
GAS=8.3144 
C      W=40. 539/4000.0 
C      W=0. 040441 
CON=GAS*TC/W 
ST=0.0 
DO  10  1=1, L 

TR=DABS ( ( T ( I ) -TO ) /TO ) 

RS= ( C ( I ) /DABS ( DSQRT ( CON ) ) ) - ( Al *TR*  *BETA+A2  *TR*  * ( BETA 
C+DELTA) +A3*TR** ( 1 . 0-ALPHA+BETA) + (A4*TR) + (A5*TR**2 . 0 ) 
C+(A6*TR**3.0) ) 
C         WRITE ( * , * )  ' RS= ' , RS 
RS=RS**2 
ST=ST+RS 
10    CONTINUE 
RETURN 
END 
C 
C 

c 

SUBROUTINE  REFLEX  (SSRSD,B,L, T,C,W) 
IMPLICIT  REAL*8  (A-H,0-Z) 
REAL* 8  DSQRT 

DIMENSION  B(20,20),C(4000) ,T(4000) ,SSRSD(500) ,A(10) 
DO  10  J=l,10 
CUM=0 . 0 
DO  50  1=1,10 
CUM=CUM+B(J,I) 
50    CONTINUE 

B(J,12)=CUM/10.0 

B(J,13)=B(J,12)+(B(J,12)-B(J,11) ) 
10   CONTINUE 

DO  200  1=12,13 
350    DO  300  J=l,10 

A(J)=B(J,I) 
300    CONTINUE 
A1=A(1) 
A2=A(2) 
A3=A(3) 
A4=A(4) 
A5=A(5) 
A6=A(6) 
TC=A(7) 
ALPHA=A(8) 
BETA=A(9) 
DELTA=A(10) 
C       IF  (I.EQ. 13.AND.TC.LT. 292. 15)  THEN 
C        WRITE (*,*)  'LOWER  BOUNDARY  LOWER  BOUNDARY' 
C         TC=292.15 

C         FACTOR=(TC-B(7,ll) ) / (B( 7 , 13 ) -B( 7 , 11 ) ) 
C         DO  600  J=l,10 
C  B(J, 13)=FACTOR*(B(J,13)-B(J, 11) )+B(J, 11) 
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C  600     CONTINUE 

C         GOTO  350 

C         ELSE  IF  (I.EQ.13.AND.TC.GT.294.15)  THEN 

C        WRITE (*,*)  'UPPER  BOUNDARY  UPPER  BOUNDARY' 

C         TC=294.15 

C  FACTOR=(TC-B(7/ll) ) / (B( 7 , 13 ) -B( 7 , 11 ) ) 

C         DO  700  J=l,10 

C  B(J,13)=FACTOR*(B(J,13)-B(J,ll) )+B(J,ll) 

C  700     CONTINUE 

C         GOTO  350 

C        END  IF 

CALL  RESIDUAL (Al,A2, A3, A4,A5,A6,TC, ALPHA, BETA 
C  , DELTA, C,T,L, ST, W) 
SSRSD(I)=ST 
C        WRITE (*,*)  'SSRSD' ,1, '=' ,SSRSD(I) 
200   CONTINUE 
RETURN 
END 
C 
C 
C 

SUBROUTINE  SHELL (B, SSRSD, WORST) 
IMPLICIT  REAL*8  (B,D-H,0-Z) 
DIMENSION  SSRSD(500) ,B(20,20) 
INTEGER  A,C 
WORST=SSRSD(ll) 
C       WRITE ( * , * )  ' SHELL ' 
J=l 
N=ll 
M=N 
90    IF  (M.GE.l)  THEN 
J=l 
M=M/2 
K=N-M 
50     IF  (J.EQ.l)  THEN 
J=0 

DO  60  I=1,K 
L=I+M 

IF  ( SSRSD ( I ) . GT . SSRSD ( L ) )  THEN 
V=SSRSD(I) 
SSRSD(I)=SSRSD(L) 
SSRSD (L)=V 
DO  70  C=l,10 
W=B ( C , I ) 
B(C,I)=B(C,L) 
B ( C , L ) =W 
70        CONTINUE 
J=l 
END  IF 
60      CONTINUE 
GOTO  50 
END  IF 
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C 


GOTO  90 
END  IF 
DO  900  1=1,11 
C        WRITE (*,*)  SSRSD(l) 
C   900  CONTINUE 
RETURN 
END 


C 
C 
C 


10 
30 


SUBROUTINE  EXPAND ( B , SSRSD , C , T , L , W ) 

IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  B(20,20) ,C(4000) ,T(4000) ,SSRSD(500) ,A(10) 

DO  10  1-1,10 

B(I,14)=B(I,13)+(B(I,12)-B(I,11)) 
CONTINUE 
DO  30  K= 1 . 1 0 


DO  30  K=l,10 

A(K)=B(K,14) 
CONTINUE 
A1=A(1) 
A2=A(2) 
A3=A(3) 
A4=A(4) 
A5=A(5) 
A6=A(6) 
TC=A(7) 
ALPHA=A(8) 


C 

C  500 


C 
IF 


ALiFHA=A(0) 

BETA=A(9) 

DELTA=A(10) 

WRITE ( * , 500 )  Al , A2 , A3 , A4 , A5 , A6 , TC, ALPHA, BETA, DELTA 

FORMAT  (1X,7F8.2,3F4.2) 
CALL  RESIDUAL (A1,A2, A3, A4,A5,A6,TC, ALPHA, BETA 
, DELTA, C,T,L, ST, W) 
SSRSD ( 14 )=ST 


20 


50 


, DELTA, C,T,L, ST, W) 
SSRSD ( 14 )=ST 

WRITE ( * , * )  ' SSRSD ■ , I , ' = • , SSRSD ( I ) 
:F  (SSRSD(14) .LT.SSRSD(13) )  THEN 
SSRSD(11)=SSRSD(14) 
DO  20  J=l,10 

B(J,11)=B(J,14) 
CONTINUE 

JLSE  IF  (SSRSD(14) .GT.SSRSD( 13 ) )  THEN 
SSRSD(11)=SSRSD(13) 
DO  50  J-1,10 

B(J,11)=B(J,13) 
CONTINUE 


C 
C 
C 
C 


END  IF 
RETURN 
END 
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SUBROUTINE  CNTR(B,SSRSD, C,T, L,W) 

IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  B(20,20) ,C(4000) ,T(4000) ,SSRSD(500) ,A(10) 

DO  10  1=1,10 

B(I,14)=B(I,12)+0.5*(B(I,12)-B(I,11) ) 
10    CONTINUE 

DO  30  K=l,10 
A(K)=B(K,14) 
30     CONTINUE 
A1=A(1) 
A2=A(2) 
A3=A(3) 
A4=A(4) 
A5=A(5) 
A6=A(6) 
TC=A(7) 
ALPHA=A(8) 
BETA=A(9) 
DELTA=A(10) 
C       WRITE(*,500)  Al ,  A2  , A3, A4,A5,A6,TC, ALPHA, BETA, DELTA 
C  500    FORMAT  ( IX, 7F8 . 2 , 3F4 .2 ) 

CALL  RESIDUAL (A1,A2, A3, A4,A5,A6,TC, ALPHA, BETA 
C  , DELTA, C,T,L, ST, W) 
SSRSD(14)=ST 
SSRSD(11)=SSRSD(14) 
DO  20  J=l,10 

B(J,11)=B(J,14) 
20    CONTINUE 
RETURN 
END 


C 
C 
C 


SUBROUTINE  CNTW(B, SSRSD,C,T,L,W) 

IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  B(20,20) ,C(4000) ,T(4000) ,SSRSD(500) ,A(10) 

DO  10  1=1,10 

B(I,14)=B(I,12)-0.5*(B(I,12)-B(I,11)) 
10    CONTINUE 

DO  30  K=l,10 
A(K)=B(K,14) 
30     CONTINUE 
A1=A(1) 
A2=A(2) 
A3=A(3) 
A4=A(4) 
A5=A(5) 
A6=A(6) 
TC=A ( 7 ) 
ALPHA=A ( 8 ) 
BETA=A(9) 
DELTA=A(10) 


C       WRITE(*,500)  A 1, A2, A3, A4,A5,A6,TC, ALPHA, BETA, DELTA 
C  500    FORMAT  ( IX, 7F8 . 2 , 3F4 .2 ) 

CALL  RESIDUAL (A1,A2, A3, A4,A5,A6,TC, ALPHA, BETA 
C  , DELTA, C,T,L, ST, W) 
SSRSD(14)=ST 
SSRSD(11)=SSRSD(14) 
DO  20  J=l,10 

B(J,11)=B(J,14) 
20    CONTINUE 
RETURN 
END 

Program  II.  Speed  of  Sound  Calculations 

C234567 

IMPLICIT  REAL*8  (A-H,0-Z) 

REAL* 8  DABS,DSQRT 

DIMENSION  C(4000),T(4000) ,A(12) 

OPEN  (60,FILE=,COOR.OUT' , STATUS= ' OLD ' ) 

DO  70  1=1,10 
READ  (60,*)  A(I) 
C        WRITE (*,*)  A(I) 
70    CONTINUE 

CLOSE (60) 

A1=A(1) 

A2=A(2) 

A3=A(3) 

A4=A(4) 

A5=A(5) 

A6=A(6) 

TC=A(7) 

ALPHA=A(8) 

BETA=A ( 9 ) 

DELTA=A(10) 

OPEN  (10, FILE= ' DATA . RAW ' , STATUS= ' OLD ' ) 

OPEN  (40, FILE= ' TC . DAT ' , STATUS= ' NEW ' ) 

L=l 
30    READ(10,*,END=20)  T(L) 
c       WRITE(*,*)  T(L) 

L=L+1 

GOTO  30 
20    CLOSE (10) 

L=L-1 

GAS=8.3144 
C       W=0. 040441 

W=0. 03007 

CON=GAS*TC/W 

DO  50  1=1, L 

TR=DABS( (T(I)-TC)/TC) 
C         WRITE (*,*)  'TR=',TR 

C(I)=DSQRT(CON)*(Al*TR**BETA+A2*TR**(BETA 
C+DELTA)+A3*TR**(1.0-ALPHA+BETA)+A4*TR+A5*TR**2.0 
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C+A6*TR**3.0) 
WRITE (40,*) 
:        WRITE (*,*) 
50    CONTINUE 
CLOSE (40) 
STOP 
END 


T(I) 
T(I) 


C(I) 
C(I) 


Program  III.  Crossing-Point  Searching 


C234567 

program  min 

CHARACTER* 20  U,V 

DOUBLE  PRECISION 

DOUBLE 

DOUBLE 

DOUBLE 

DOUBLE 

DOUBLE 

DOUBLE 

DOUBLE 


Q,T,S,R,E,F 

DEL ,DSQRT, DABS, ZERO, TWO, DELE 

P2L,P3L,TCL,TL,TRL 

P2R,P3R,TCR,TR,TRR,DELF 

A1L,A2L,A3L,A4L,A5L,A6L 

ALPHAL , BETAL , DELTAL 

A1R,A2R,A3R,A4R,A5R,A6R 

ALPHAR , BETAR , DELTAR 

B(20) 


PRECISION 
PRECISION 
PRECISION 
PRECISION 
PRECISION 
PRECISION 
PRECISION 
DIMENSION  A(20) , 
ZERO=0.0D00 
TWO=20.0D-1 

WRITE (*,*)  'RIGHT  SIDE  FILE  NAME 
READ(*,50)  U 
50    FORMAT  (A20) 

OPEN  (20, FILE=U , STATUS= ' OLD ' ) 
DO  40  1=1,10 
READ  (20,*)  A(I) 
WRITE (*,*)  A(I) 
40    CONTINUE 
CLOSE (20) 

WRITE (*,*)  'LEFT  SIDE  FILE  NAME' 
READ(*,50)  V 

OPEN  ( 3  0 , FILE=V , STATUS= ' OLD ' ) 
DO  60  1=1,10 
READ  (30,*)  B(I) 
WRITE (*,*)  B(I) 
60    CONTINUE 
CLOSE (30) 
A1R=A(1) 
A2R=A(2) 
A3R=A(3) 
A4R=A(4) 
A5R=A(5) 
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A6R=A(6) 
TCR=A(7) 
ALPHAR=A ( 8 ) 
BETAR=A(9) 
DELTAR=A(10) 
A1L=B(1) 
A2L=B(2) 
A3L=B(3) 
A4L=B(4) 
A5L=B(5) 
A6L=B ( 6 ) 
TCL=B ( 7 ) 
ALPHAL=B ( 8 ) 
BETAL=B ( 9 ) 
DELTAL=B(10) 

WRITE (*,*)  'INPUT  TR  &  TL  IN  KELVIN' 
READ(*,*)  TR,  TL 
T=(TR+TL)/TWO 
P2L  =  BETAL  +  DELTAL 
P3L  -  1.0-ALPHAL+BETAL 
P2R  =  BETAR  +  DELTAR 
P3R  =  1.0-ALPHAR+BETAR 
10    TRL  =  DABS ( ( T-TCL ) /TCL ) 
TRR  =  DABS( (T-TCR)/TCR) 
:      write  (*,*)  trl,trr,p21,p31,p2r,p3r 

S= ( DSQRT ( TCL ) ) * (A1L*TRL**BETAL+A2L*TRL**P2L 
C+A3L*TRL*  *P3L+A4L*TRL+A5L*TRL*  *2+A6L*TRL*  *  3 ) 
R= ( DSQRT ( TCR ) ) * ( AlR*TRR*  *BETAR+A2R*TRR*  *P2R 
C+A3R*TRR**P3R+A4R*TRR+A5R*TRR**2+A6R*TRR**3) 
DEL=R-S 
IF  (DABS (DEL) .LE.1.0D-16)  THEN 

GOTO  100 
!       ELSE  IF  (R.GT.S)  THEN 

ELSE  IF  ( DEL. GT. ZERO)  THEN 

TR=T 

T=(TR+TL)/TWO 

DELE=T-E 

IF  ( DELE. EQ. ZERO)  THEN 
GOTO  90 

END  IF 

E=T 

WRITE  (*,*)  DEL,T,TR 

GOTO  10 
I      ELSE  IF  (R.LT.S)  THEN 

ELSE  IF  ( DEL. LT. ZERO)  THEN 

TL=T 

T=(TR+TL)/TWO 

DELF=T-F 

IF  (DELF.EQ.ZERO)  THEN 
GOTO  95 

END  IF 

F=T 
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WRITE  (*,*)  DEL,T,TL 
GOTO  10 
END  IF 
90   WRITE  (*,*)  '****', T,E 

GOTO  97 
95   WRITE  (*,*)  '****', T,F 

GOTO  97 
97   WRITE  (*,*)  * — >UNABLE  TO  REFINE  FURTHER 
100   WRITE  (*,*)  T,T-273.15,DEL 
R=R*(SQRT(8. 3144/0. 03007) ) 
S=S*(SQRT(8. 3144/0. 03007) ) 
WRITE  (*,*)  R,S 
STOP 
END. 


APPENDIX  D 
LABORATORY  STANDARD  PRESSURE  GAUGE 


Introduction 

Pressure  gauges  are  essential  instruments  of  science  and 
industry  and  are  available  in  a  variety  of  types.  Most  of  the 
common  types  are  favored  for  possessing  one  or  more  of  the 
following  desirable  features:  low  cost,  ease  of  use,  rugged 
construction,  suitability  for  remote  sensing.  The  attainment 
of  these  features  is  generally  facilitated  by  the  use  of 
gauges  which  measure  pressure  indirectly.  Such  gauges  require 
calibration  and  are,  therefore,  not  suitable  for  use  as 
primary  standards.  Unless  their  calibrations  are  checked 
frequently,  their  readings  are  always  somewhat  suspect. 

While  for  most  purposes  the  common  gauge  types  give 
sufficiently  reliable  readings,  there  are  occasions  which  call 
for  absolute  pressure  measurements  made  by  gauges  which 
measure  directly  the  average  normal  force  acting  over  a  known 
area.  The  pressure  range  from  about  10"3  to  10"5  pascal  is 
served  adequately  at  the  low  end  by  the  McLeod  gauge  if 
suitably  operated  to  avoid  the  errors  associated  with  the 
mercury  pumping  effect142  and  at  the  upper  end  by  the  U-tube 
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148 
manometer.  Higher  pressures  require  the  use  of  pressure  gauges 
of  the  dead  weight  variety143. 

This  paper  reports  on  the  construction  and  use  of  a 
simple  dead  weight  gauge  for  absolute  pressure  measurements 
and  a  variation  of  that  gauge  which  works  well  but  requires 
calibration.  Figure  D.l  illustrates  the  basic  features  of  a 
conventional  dead  weight  pressure  gauge.  The  pressure  to  be 
measured  is  transmitted  to  the  floating  piston  through  a  null 
detector  such  as  a  mercury  U-tube  with  adjustable  electric 
contacts,  for  example.  The  null  detector  is  used  to  isolate 
the  system  gases  from  the  lubricating  oil  which  surrounds  the 
piston.  This  oil  is  continually  lost  by  viscous  flow  and  must 


Figure  D.l.  Conventional  dead  weight  pressure  gauge:  A  = 
piston;  B  =  cylinder;  C  =  steel  U  tube;  D  =  oil  injector;  N, 
N'  =  indicator  contact  needles. 
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Fxgure  D.2.  The  present  deadweight  pressure  gauge. 
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be  replenished  by  means  of  an  injection  pump.  The  upward  force 
on  the  piston  is  balanced  by  the  sum  of  the  external  pressure, 
gravity  acting  on  the  piston  and  an  additional  masses  added  as 
required  to  keep  the  piston  afloat.  Although  this  type  of 
pressure  gauge  is  capable  of  high  accuracy,  its  use  is 
tedious,  and  commercial  units  are  expensive,  being  beyond  the 
means  of  many  small  laboratories.  The  present  gauges  were 
developed  to  avoid  these  problems. 

Dead  Weight  Pressure  Gauge 

The  present  gauge  is  shown  in  figure  D.2.  It  differs 
principally  from  the  conventional  designs  in  the  following 
respects: 

1.  The  piston  is  inverted  being  forced  vertically 
downward  by  the  system  pressure. 

2.  The  loading  weight  is  generated  by  standard  masses 
suspended  from  accurately  spaced  hooks  on  a  pivoted  balance 
beam. 

3.  The  load  is  transferred  to  the  piston  by  point  contact 
with  a  hardened  precision  ball  accurately  located  on  the  beam 
with  respect  to  the  hooks. 

4.  A  thin  flexible  membrane  across  the  top  of  the 
cylinder  acts  as  a  piston  follower  thereby  separating  the 
piston  from  the  system  gases  and  eliminating  the  need  for  a 
null  detector  or  oil  injection  pump. 
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Figure  D.3  shows  the  arrangement  of  the  gauge  elements. 
System  pressure  is  delivered  to  the  top  surface  of  the 
follower  and  transmitted  through  it  to  the  top  of  the  piston 
which  slides  freely  in  the  mating  cylinder.  The  resultant 
force  is  countered  by  the  upward  force  of  the  loaded  balance 
beam.  The  travel  of  the  piston  is  limited  by  a  stop  pin 
extending  into  an  oversize  hole  in  the  beam.  Total  vertical 
displacement  of  the  piston  and  angle  between  the  beam  and  true 
horizontal  influence  the  net  downward  force.  For  this  reason 
the  balance  point  is  always  taken  as  the  median  position,  for 
which  the  beam  is  horizontal  and  the  relationship  P=F  /A 
holds,  where  Ff  is  the  upward  force  transmitted  by  the  piston 
to  the  follower.  This  gauge  was  built  for  the  specific  purpose 
of  measuring  pressures  and  standardizing  secondary  gauges  used 
in  an  ongoing  research  effort  on  sonic  speeds  in  gases  at 
pressures  up  to  about  103  pascal  (-  150  psi.)  The  stresses 
generated  on  the  gauge  elements  at  these  pressures  are 
relatively  small  so  ordinary  materials  could  be  used.  The 
piston,  cylinder  and  gauge  body  were  machined  from  type  303 
stainless  steel.  For  use  at  high  pressures  these  elements 
should  be  made  from  heat  treatable  alloys  and  hardened  in  the 
conventional  manner.  Many  materials  are  available  for  use  as 
follower.  A  choice  may  be  made  principally  on  the  basis  of 
compatibility  with  the  fluids  to  be  measured  and  flexibility. 
Low  permeability  is  desirable,  of  course,  but  since  the 
follower  is  supported  by  dense  solid  surfaces,  loss  of  gas 
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through  it  is  generally  negligible.  Some  suitable  follower 
materials  are  thin  elastomeric  sheets,  plastic  films  and  metal 
foils  if  not  stressed  beyond  their  elastic  limits.  The  balance 
beam  was  machined  from  brass  bar  stock.  The  crucial  features 
on  it  are  the  location  of  the  contact  ball  and  the  weight 
suspension  hook  pins  with  respect  to  the  knife  edge.  These 
features  were  located  on  2.000  inch  centers  with  maximum 
uncertainty  of  0.0005  inch.  The  contact  ball  sits  in  a  60° 
tapered  center  drill  hole  on  top  of  the  beam,  and  the  hardened 
0.1250  inch  hook  pins  are  pressed  into  mating  holes  located 
along  the  lower  edge  of  the  beam.  The  S  hooks  fit  over  the 
pins  and  pivot  freely  in  a  groove  machined  into  the  bottom  of 
the  beam.  The  pins  are  numbered  consecutively  one  to  ten 
beginning  with  the  one  nearest  the  fulcrum.  The  beam  is 
pivoted  on  a  knife  edge.  Conventional  balances  have  quartz  or 
other  hard  crystalline  knife  edges  and  bearing  pads.  Because 
such  components  are  expensive  and  not  readily  available, 
especially  for  supporting  heavy  loads,  it  was  decided  to  make 
use  of  a  common  low  cost  material  which  seemed  well  suited  for 
the  purpose.  Cemented  carbide  inserts  produced  by  powder 
metallurgical  techniques  for  use  as  cutting  tools  in  machining 
operations  are  available  in  many  standard  compositions  and 
geometries.  These  are  harder  than  quartz  and  are  manufactured 
to  exacting  dimensions  with  close  tolerances.  We  used  a 
standard  insert  0.500  X  0.500  X  0.125  inch  which  was  set  into 
a  seat  cut  at  a  45°  angle  with  respect  to  the  bottom  of  the 
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beam.  The  insert  has  eight  useable  90°  knife  edges.  A  second 
insert  was  used  as  the  bearing  pad.  It  fits  in  a  seat  machined 
in  a  support  member  and  securely  fastened  to  the  gauge  frame. 
When  not  in  use  the  beam  may  be  elevated  to  lift  the  knife 
edge  off  of  the  bearing  pad  to  prolong  the  life  of  the  knife 
edge  by  supporting  it  on  the  hardened  pins. 

The  cylinder  was  machined  by  drilling  and  reaming  and 
finished  by  ballizing.  The  piston  was  machined  on  a  precision 
lathe  and  polished  with  rouge  paper.  The  diameter  was  measured 
with  a  Pratt  and  Whitney  supermicrometer .  Ten  readings  taken 
at  different  positions  around  the  piston  gave  the  average 
diameter  0.37383  inch  with  a  standard  deviation  of  0.00001 
inch. 

Because  in  operation  the  contact  ball  moves  on  a  circular 
arc  and  the  piston  is  constrained  to  move  only  vertically, 
there  is  some  relative  motion  between  them  resulting  in  some 
unwanted  sliding  friction.  To  avoid  this  problem  a  small 
thrust  bearing  was  used  between  the  lower  surface  of  the 
piston  and  a  hardened  chip  to  convert  the  sliding  motion  to 
the  less  restricted  rolling  motion. 

The  gauge  is  mounted  on  a  rigid  beam  which  is  securely 
bolted  to  the  laboratory  floor.  The  mounting  plate  is 
adjustable  to  allow  the  bearing  pad  to  be  precisely  leveled. 
The  vertical  component  of  upward  force  exerted  by  the  contact 
ball  on  the  gauge  piston  is  equal  to 

Fr  =  g  cos  0  S  n^  (D«l) 
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where  g  =  local  acceleration  of  gravity, 

m.  =  added  load  mass 

ni  =  position  of  load  mass  and 

0   =  the  angle  between  the  top  of  the  beam  and  the 
horizontal. 
In  order  to  use  the  simpler  relationship 

FT  =  g  S  n.m.  (D.2) 

it  is  necessary  to  adjust  the  position  of  the  gauge  body 
relative  to  the  plane  of  the  bearing  pad  so  that  the  balance 
point  occurs  when  the  beam  is  precisely  horizontal  (i.e.  cos 
0  =  1.)  To  accomplish  this  the  gauge  body  is  mounted  in  a 
grooved  way  and  adjusted  vertically  by  a  fine  lead  screw.  A 
precision  spacer  is  positioned  between  the  lower  end  of  the 
cylinder  and  the  flange  on  the  piston  during  the  adjustments. 
This  spacer  assures  that  the  piston  can  extend  into  the 
cylinder  only  to  the  median  position.  The  beam  is  loaded 
sufficiently  to  hold  the  piston  against  the  spacer  and  the 
gauge  body  adjusted  vertically  by  means  of  the  lead  screw  so 
that  a  precision  level  shows  the  beam  to  be  horizontal. 
Following  the  adjustments  the  gauge  body  is  clamped  securely 
in  place. 

While  the  piston  is  held  in  the  median  position  by  use  of 
the  spacer,  the  balance  scale  and  indicator  are  set  to  "zero". 
Although  a  conventional  scale  and  pointer  combination  can  be 
used  with  this  apparatus,  greater  sensitivity  is  easily 
achieved  by  using  an  optical  indicator  consisting  of  a  He-Ne 
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laser  beam  reflected  off  a  mirror  attached  to  the  balance  beam 
onto  a  suitable  scale.  We  have  used  a  sheet  of  precision  graph 
paper  taped  to  a  wall  of  the  laboratory.  The  path  length 
between  the  mirror  and  the  scale  is  about  3  meters . 

The  unloaded  beam  is  first  balanced  by  using  the 
adjusting  weights  added  as  necessary  to  restore  the  balance  of 
forces  on  the  piston.  To  reduce  the  effects  of  friction,  which 
can  lead  to  weighing  errors,  the  piston  is  coated  with  oil  and 
provided  with  a  horizontal  lever  arm  which  is  used  to 
oscillate  it  back  and  forth  through  a  small  angle  while  it 
settles  into  its  vertical  eguilibrium  position.  Freedom  from 
sticking  is  evident  when  the  balance  point  is  constant  on 
repeated  tries  including  those  with  both  positive  and  negative 
initial  offsets. 
Test  Results 

The  absolute  dead  weight  pressure  gauge  was  tested  for 
accuracy  by  using  it  to  determine  the  pressure  of  gaseous 
CF2C12  (freon-12)  in  eguilibrium  with  its  liquid  phase  at 
constant  temperature.  This  system  has  been  thoroughly  studied 
by  Gilkey  et.  al.(144'  145)  using  more  elaborate  methods.  The 
freon  container  was  placed  in  an  adjustable  subambient 
constant  temperature  bath  and  connected  to  the  gauge  inlet. 
Air  was  purged  from  the  gauge  and  lines  by  vacuum  pumping  and 
repeated  flushing  with  small  aliquots  of  freon.  After 
equilibrium  was  reached  the  beam  was  balanced  and  the 
requisite  masses  and  their  positions  noted.  The  atmospheric 
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pressure  was  determined  using  a  mercury  barometer  noting  the 
usual  corrections  for  thermal  expansions  of  mercury  and  brass. 
The  measured  pressure  was  calibrated  from  the  relationship 

P  "  pata,  +  (9  s  n.m.)/A  (D.3) 

where  A  =  nr2  =  cross  sectional  area  of  the  piston 
The  local  acceleration  of  gravity  must  be  known  or  measured. 
In  this  laboratory  it  is  979.284  cm/sec.2  and  the  resulting 
gauge  constant  is  0.0200584  psig/g. 

Table  D.l  gives  the  results  of  these  tests  at  two 
temperatures,  along  with  the  corresponding  vapor  pressures 
which  may  be  inferred  from  the  precise  measurements  of  Gilkey 
et.  al.<144/  145)  The  close  agreement  shows  that  the  present  gauge 
works  well.  These  results  are  not  intended  to  serve  as  a 
calibration  of  the  gauge  but  rather  as  a  check  of  its 
performance. 

The  f reon  used  was  manufactured  for  refrigeration  service 
and  was  of  unknown  purity.  The  consistency  of  the  gauge 
readings  indicates  that  the  absolute  pressure  gauge  is 
accurate  to  within  the  limits  of  the  smallest  mass  difference 
which  may  be  reproducibility  discerned.  Using  a  30  ^m  polyimid 
film  follower  (DuPont  Kapton)  the  loading  mass  could  be 
measured  easily  to  within  ±  0.5  g,  which  corresponds  to  the 
limits  of  uncertainty  stated  in  the  table.  The  gauge  is  a 
little  more  responsive  if  a  15  ^m  latex  follower  is  used  (ca. 
±  0.1  g)  but  the  lubricating  oils  used  thus  far  permeate  or 
react  with  the  latex  causing  operational  difficulties.  Thin 
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films  of  synthetic  elastometers  such  as  neoprene,  viton, 
silicone  rubbers  etc.  should  work  well,  but  have  not  thus  far 
been  available  for  testing. 

The  tests  included  comparison  of  the  readings  of  a 
precision  Bourdon  gauge  connected  to  the  pressure  line.  Table 
D.l.  shows  that  this  particular  gauge  is  in  error  or  requires 
calibration. 


Deadweight  Gauge  with  Bellows  Sensor 
The  basic  pressure  gauge  herein  described  was  designed  to 
accommodate  interchangeable  sensing  elements.  It  can,  for 
Table  D.l.  Experimental  Test  of  Dead  Weight  Pressure 
Gauge . 


Temperature,  K 

Vapor  Pressure  R-12,  psia 

(in  parenthesis  is  %  deviation 

from  literature  value) 

Gilkey 
(1947) 

This  work 

Bourdon 
Gauge 

271.92 

43.04 

42.97±0.01 
(0.16%) 

42.47 

(1.3%) 

287.45 

69.87 

69.66 
(0.31%) 

68.90 
(1.4%) 

example,  be  fitted  with  piston/cylinder  combinations  of 
different  diameters  or  with  entirely  other  types  of  sensors. 
In  particular  it  was  intended  to  investigate  the  suitability 
of  a  thin  walled  metallic  bellows  for  use  as  a  balance 
detector.  An  electrof ormed  nickel  bellows  with  0.0018  inch 
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walls,  0.375  inch  o.d.,  0.25  inch  i.d.  and  24  convolutions  was 
brazed  to  fittings  which  permitted  it  to  be  utilized  as  a  re- 
entrant pressure  sensor  (see  figure  D.4.)  System  pressure  is 
applied  externally  around  the  bellows  causing  it  to  contract 
or  shorten.  A  push  rod  extending  through  the  bellows  is  acted 
upon  by  the  contact  ball  as  before.  The  beam  is  loaded  to 
counteract  the  downward  force  due  to  the  system  pressure.  The 
bellows  are  very  flexible  axially  but  stiff  in  the  radial 
direction.  The  downward  force  is  then  equal  to  the  magnitude 
of  the  system  pressure  times  the  active  cross  sectional  area, 

Aeff  °r 

P  -  VAeff  "  p.tm  +  g2mini/Aeff.  (D.4) 

Ae£f  is  given  approximately  by  the  bellows  manufacturer  as 
0.0723  inch2,  but  lack  of  accuracy  in  this  value  and 
uncertainty  in  the  effects  of  the  brazing  operations  on  it 
require  that  it  be  determined  by  calibration.  Freon-12  was 
used  for  calibrating  the  bellows.  The  weight  required  to 
balance  the  system  pressure  was  measured  versus  temperature  of 
the  vapor-liquid  equilibrium.  If  the  bellows  does  behave  as  a 
piston  of  fixed  active  cross  section,  a  plot  of  system 
pressure  versus  the  mass  required  to  balance  the  beam  should 
be  a  straight  line  with  slope  of  g/Aeff.  Figure  D.5  shows  a 
plot  of  the  calibration  data.  The  slope  was  determined  by  a 
linear  regression  and  is  equal  to  0.03000273  psig/g.  The 
correlation  coefficient  of  the  linear  fit  is  0.99995  which 
shows  the  gauge  to  perform  in  accordance  with  equation  D.4 
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Figure  D.4.  An  electrof ormed  nickel  bellow  pressure  sensor, 
over  the  range  of  pressures  tested.   The  effective  area 
inferred  by  the  calibration  is  0.0734  inch2. 

The  disadvantage  of  this  gauge  relative  to  the  piston 
type  is  the  need  for  calibration.  The  advantage  is  the 
intrinsic  greater  sensitivity  available  with  the  bellows.  The 
present  sensor  responds  to  mass  differences  of  0.05  gram  which 
corresponds  to  pressure  differences  of  0.0015  psi. 

Conclusion 

This  work  shows  that  even  small  laboratories  can  have 
reliable  pressure  standards  at  relatively  low  cost.  The  gauges 
described  here  were  constructed  from  readily  available 
materials  at  a  total  cost  of  less  than  $300  (excluding  the 
laser.)  Commercial  dead  weight  pressure  gauges  cost  several 
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thousand  dollars  and  are  generally  beyond  the  means  of  many  of 
the  laboratories  which  could  benefit  from  their  use.  Perhaps 
gauges  based  on  the  technology  presented  here  will  bridge  the 
present  gap  between  conventional  secondary  gauges  and  high 
cost  absolute  ones.  They  are  at  least  working  well  in  this 
laboratory  by  increasing  the  confidence  in  the  determination 
of  pressure  sensitive  acoustic  and  thermodynamic  properties. 
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APPENDIX  E 

COMPUTER  PROGRAMS 

Since  most  of  the  programs  are  lengthy  only  programs 

VIKING  and  MAX3  are  showed  in  full.  The  others  are  listed  only 

main  programs.  Most  of  these  programs  were  written  in  fortran 

77  by  Dr.  McGill.  The  remaining  by  the  author. 


PROGRAM  VIKING 

$ large 

INTERFACE  TO  SUBROUTINE  TIME  (N,STR) 
CHARACTER* 10  STR  [NEAR, REFERENCE] 


INTEGER* 2  N  [VALUE] 
END 


REAL*4  DX(2000) ,DY(2000) , XMAX , XMIN , YMAX , YMIN 
CHARACTER* 1  KOL 
CHARACTER* 2  SWEEP 
CHARACTER* 10  DNAME 
CHARACTER* 7  DTEMP 
CHARACTER* 10  TSTR 
INTEGER*2  HR,MIN,SEC 
INTEGER  READS, T,TC ,10,100, POINTS 
SWEEP='30' 
WIDTH  =8.0 
READS  =  124 
KOUNT  =  0 
STP  =  0.0025 
CALL  KEITHA(TZERO) 

OPEN ( 6  0 , File= ' c : DATA . BUF  * , STATUS= ' OLD ' ) 
READ (60,*)  VOLT,IG 
CALL  GAINSET(IG) 
CLOSE (60) 
C*********************************************************** 

C       VOLTE  =  (25.0-TZERO)*0.01  +  VOLT 
C  TO  BE  USED  IN  NORMAL  SCAN 
C********************************* ************************** 

VOLTE  =  0.0175 
CLOSE (10) 
DTEMP='c:DATA. ■ 

OPEN (70, File= ' c : NAME . NUM ' , STATUS= ' OLD ' ) 
READ(70,301)  DNAME 
301      FORMAT ( IX, A10) 
CLOSE  (70) 

READ(DNAME, * (7X,I3) ' )  NUM 
NUM=NUM+1 
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OPEN (70, File= ' c : NAME . NUM • , STATUS= ' NEW ' ) 
WRITE (70, 302)  DTEMP,NUM 
WRITE(*,302)  DTEMP,NUM 
302      F0RMAT(1X,A7,I3) 
CLOSE  (70) 

OPEN (20, FILE=DNAME , STATUS= ' NEW ' ) 
OPEN(30,File=,c:Peak.buf , STATUS= ■ OLD ' ) 
READ(30,*)  F,TI 
TI=TI+273.15 
CALL  KEITHA(TF) 
TF=TF+273.15 
F=SQRT(TF/TI)*F 
CLOSE (30) 
IXT=3 
c  IXT=1 

NA  =  21/IXT 
IDELY=2*NA 

Q*  ********************************************************  * 

C  M=1000/NA/IXT-IDELY/NA 

C  USED  ONLY  FOR  CENTERING 

Q********************************************************** 

M=61 
C  USED  AFTER  CENTERING 

Q********************************************************** 

XMIN=1.0 
XMAX=FLOAT(M) 
1001     CALL  TIME(10,TSTR) 

READ(TSTR, ' (12) ',ERR=1001)  HR 
READ(TSTR, ' (3X,I2) ' ,ERR=1001)  MIN 
READ(TSTR, ■ (6X,I2) ' ,ERR=1001)  SEC 
TOO  =  HR*3600+MIN*60+SEC 

100  CALL  TEMPSET(VOLT) 

101  DO  30  1=1, READS 
500        SF=F-WIDTH/2.0 

SPF=F+WIDTH/2.0 
INDEX=1 
POINTS=M 
ISTART=1 
200        CALL  HP (SF,SPF, SWEEP) 
AMAX=0 . 0 
XMIN=200000.0 
DO  42  ID=1,IDELY 
CALL  KEITH (MV) 
42         CONTINUE 

DO  40  J=ISTART, POINTS, INDEX 
AVG=0 . 0 
DO  41  JA=1,NA 
300  CALL  KEITH (MV) 

KOUNT=KOUNT+l 
IF(KOUNT.EQ.2*NA)  THEN 
IG=IG+1 
CALL  GAINSET(IG) 
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ENDIF 

IF(K0UNT.GT.3*NA)  THEN 
KOUNT=0 
GOTO  500 
ENDIF 

IF(MV. GT.20000.OR.MV.LT. 100)  GOTO  300 
KOUNT=0 

AVG=AVG+ FLOAT ( MV ) 
4 1         CONTINUE 

AMP=AVG/ FLOAT ( NA ) 
DX ( J ) =FLOAT ( J ) 
DY(J)=AMP 
IF(AMP.GT.AMAX)  THEN 
AMAX=AMP 
JMAX=J 
ENDIF 

IF(AMP.LT.XMIN)  XMIN=AMP 
40         CONTINUE 

IF ( INDEX. EQ.l)  THEN 

NPLOT=l 
ELSE 

NPLOT=0 
ENDIF 
YMAX=AMAX 

Q*  *****************************************************  * 

C  CALL  PLOT(M,DX,DY,XMAX,XMIN,YMAX,YMIN,NPLOT) 

C  ONLY  USED  FOR  CENTERING 

C*********** ******  *********  ***************************** 

IF (INDEX. EQ.l)  THEN 

INDEX=-INDEX 

TEMP=SF 

SF=SPF 

SPF=TEMP 

ISTART=POINTS 

POINTS=l 

JFMAX=JMAX 

AVGA=AMAX 

CALL  KEITHA(TI) 

IF(TI.LT.25.0)  GOTO  169 

CALL  PRES(VI) 

GOTO  200 
ELSE 

CALL  KEITHA(TF) 

CALL  PRES(VF) 

JRMAX=JMAX 

FO=FLOAT ( JFMAX+ JRMAX ) / 2 . 0 

CALL  MAX(F, WIDTH, FO,M) 

F=FO 

AMP= ( AMAX+AVGA ) / 2 . 0 

IF ( JFMAX . EQ . 1 . OR . JFMAX . EQ . M . OR . JRMAX . EQ . 1 . OR . 
C JRMAX. EQ.M)  THEN 

OPEN(30,File='c:Peak.buf ' , STATUS= ' OLD ■ ) 
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READ (30,*)  F,TI 
CLOSE (30) 
TI=TI+273.15 
TF=TF+27  3.15 
F=SQRT(TF/TI)*F 
GOTO  500 
ENDIF 
1003  CALL  TIME(10,TSTR) 

READ(TSTR, ' (12) ' ,ERR=1003)  HR 
READ(TSTR, ' (3X,I2) ' ,ERR=1003)  MIN 
READ(TSTR, ' (6X,I2) ' ,ERR=1003)  SEC 
T  =  HR*3600+60*MIN+SEC 
T=T-T00 

IF(T.LT.O.O)  THEN 
T00=T00-24*3600 
T  =  HR*3600+60*MIN+SEC 
T=T-T00 
ENDIF 
IF(AMP.LT.7000)  THEN 
IG=IG-1 

IF(IG.GE.7)  CALL  GAINSET(IG) 
ENDIF 

IF(AMP.GT. 19000)  THEN 
IG=IG+1 

IF(IG.LE.24)  CALL  GAINSET(IG) 
ENDIF 

RATIO=F/SQRT ( ( TI+TF+273 . 15*2 . 0 ) /2 . 0 ) 
WRITE(20,*)  TI,TF,VI,VF,AMP,F,T,IG 

Q*  *********************************************************  * 

C  DO  18  IF=1,19 

C  WRITE ( * , * ) 

C  18          CONTINUE 

C  WRITE (*, 390)  SPF,SF 

C  390         FORMAT ( IX, ' FREQUENCY=  ' ,F10 . 3 , 48X,F10 . 3 ) 

C  WRITE ( * , * ) 

C  IF  (JRMAX  .GT.  JFMAX  )  M=M-1 

C  IF  (JRMAX  .LT.  JFMAX  )  M=M+1 

C  DF=0 . 0 

C  USED  ONLY  FOR  CENTERING 

Q*********************************************************** 

WRITE (*, 400 )TI,TF, VI, VF,M,F,T, JRMAX, JFMAX, RATIO 

400       FORMAT(1X,2F8.3,2F8.1,I6,F10.3,I8,2I4,F9.3) 
(2*  *********************************************************  * 

DT=TF-TI 

DF=FLOAT ( JRMAX- JFMAX ) /FLOAT ( M ) *WIDTH 
C  USED  ONLY  AFTER  CENTERING 
(2*********************************************************** 

F=F+DF*1.5 

ENDIF 
30    CONTINUE 
C        IF(ABS(DT) .GT.0.02)  READS=READS+1 
C        IF(ABS(DT) .GT.0.03)  READS=READS+1 
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C 
C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


IF(ABS(DT 
IF(ABS(DT 
IF(ABS(DT 
IF(ABS(DT 
IF(ABS(DT 
IF(ABS(DT 
IF(ABS(DT 
IF(ABS(DT 
IF(ABS(DT 
IF(ABS(DT 
IF(ABS(DT 
IF(ABS(DT 


.GT.0.04) 

.GT.0.05) 

.GT.0.06) 

.GT.0.07) 

.GT.0.08) 

.LT.0.05. 

.LT.O 

.LT.O 

•  LT.O 
.LT.O 

•  LT.O 


READS=READS+1 

READS=READS+1 

READS=READS+1 

READS=READS+1 

READS=READS+1 

AND . READS . GT . 2 

04. AND. READS. GT. 2 

03. AND. READS. GT. 2 

02. AND. READS. GT. 2 

015. AND. READS. GT. 

0  LAND.  READS.  GT.  2 


)  READS=READS-1 
)  READS =RE ADS -1 
)  READS=READS-1 
)  READS=READS-1 
2)  READS=READS-1 
)  READS=READS-1 
2)  READS=READS-1 


380 


C****** 

c 
c 

c 
c 
c 
c 
c 
c 
c 

C  USED 

C****** 
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C 
C 

c 


. LT.O. 005. AND. READS. GT. 
WRITE (*, 380)  VOLT 

FORMAT  (IX,' SET  POINT  VOLTAGE  =  ' ,F9.4) 
OPEN(30,File='c:Peak.buf ' , STATUS= ' NEW ' ) 
WRITE (30,*)  F,TF 
CLOSE (30) 

IF(TF.GT. 100.0)  GOTO  169 
IF ( VOLT. LE. VOLTE)  THEN 

VOLT=VOLT+STP 

OPEN(60,File='c:VOLT.BUF' , STATUS= ' NEW ' ) 

WRITE (60,*)  VOLT 

CLOSE (60) 

GOTO  100 
********************************************* 

ELSE 

IF ( VOLT. GT. VOLTE- 0.4)  THEN 

VOLT=VOLTE-0.5 

OPEN ( 60 ,File='c: VOLT. BUT' , STATUS= ' NEW • ) 

WRITE(60,*)  VOLT 

CLOSE (60) 

GOTO  100 

ENDIF 

GOTO  101 

IN  NORMAL  SCAN 
******************************************************* 

ENDIF 

OPEN ( 70 ,File='c: STOP. RUN' , STATUS= ' NEW ' ) 

CLOSE (70) 

CALL  TEMPSET(O.O) 

STOP 

END 


SUBROUTINE  KEITHA(T) 

common  /ibglob/  ibsta, iberr , ibcnt 

' COMMON  GROUP  1 . ' 

integer*4      cmd(10) rrd(512) ,wrt(512) 

character* 8   bname , bdname , TEMP 

1  COMMON  GROUP  2 . ■ 

character*50  flname 
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KLA  =  IBFIND  ('K195A  ') 
124  =  2**24 
116  =  2**16 

18  =  2**8 
1004     CALL  IBRD  (KLA,RD,12) 

15  =  RD(2)/I24 

16  =  (RD(2)-I5*I24)/I16 

17  =  (RD(2)-I5*I24-I6*I16)/I8 
188  =  RD(2)-I5*I24-I6*I16-I7*I8 

19  =  RD(3)/I24 

110  =  (RD(3)-I9*I24)/I16 

111  =  (RD(3)-I9*I24-I10*I16)/I8 

112  =  RD(3)-I9*I24-I10*I16-I11*I8 
WRITE(TEMP, • (8A1) ■ )char(l88) ,char(I7) ,char(l6) ,char(I5) , 

CCHAR(I12) ,CHAR(Ill) , CHAR (I 10) , CHAR (19) 
READ(TEMP, • (F8.3) ' ,ERR=1004)  R 
R0=99.99 
ALPHA=0. 0039042 
DELTA=1.5205 
c  R0=100.06 

c  ALPHA=. 0039046 

c  DELTA=1.5205 

ALDEL=ALPHA* DELTA 
RC=R/R0-1.0 

PART1=ALPHA+ALDEL/ 100.0 
PART2= ( ALDEL/100 . 0+ALPHA) **2 
PART3=4 . 0*RC*ALDEL/10000 . 0 
PART4=2 . 0*ALDEL/10000 . 0 
TEMP2=PART1-SQRT(PART2-PART3 ) 
T=TEMP2/PART4 
RETURN 
end 
C 
C 
C 

SUBROUTINE  HP (SF,SPF, SWEEP) 
common  /ibglob/  ibsta, iberr,ibcnt 
' COMMON  GROUP  1 . ' 

integer*4      cmd(10) ,rd(512) ,wrt(512) ,IP(4) 
character*l  ch(40) 
CHARACTER* 9  STARTF,STOPF 
CHARACTER* 2  SWEEP 
character* 8   bname,bdname 
' COMMON  GROUP  2 . ' 
call  ibinit (ibsta) 
IP(1)  =  1 
IP(2)  =  2**8 
IP(3)  =  2**16 
IP(4)  =  2**24 
IHP  =  IBFIND  ('hp3325  ') 
CH(1)  =  'F' 
CH(2)  =  'U' 
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CH(3)  -  ' 1' 

CH(4)  =  'A' 

CH(5)  =  'M' 

CH(6)  =  '0' 

CH(7)  =  'l' 

CH(8)  =  'V 

CH(9)  =  '0' 

CH(10)  =  'S' 

CH(ll)  =  'T' 

WRITE ( STARTF , ■ ( F9 . 3 ) 

READ (STARTF, ■ (Al) • ) 

READ (STARTF, ■ (1X,A1 

READ (STARTF, ' (2X,Al 

READ (STARTF, ' (3X,Al 

READ (STARTF, ' (4X,Al 

READ (STARTF, ' (5X,A1 

READ (STARTF, ' (6X,A1 

READ (STARTF, ' (7X,Al 

READ (STARTF, ' (8X,Al 

READ (STARTF, ' (9X,Al 

CH(22)  =  'H' 

CH(23)  =  'Z' 

CH(24)  =  'S' 

CH(25)  =  'P' 

WRITE (STOPF, ' 

READ(STOPF, ■ ( 

READ (STOPF, ' ( 

READ (STOPF, ' ( 

READ (STOPF, ■ ( 

READ (STOPF, ' ( 

READ (STOPF, ' ( 

READ (STOPF, ' ( 

READ (STOPF, ' ( 

READ (STOPF, ' ( 

READ (STOPF, ' ( 

CH(36)  =  'H' 

CH(37)  =  'Z' 

CH(38)  =  'T' 

CH(39)  =  '!■ 

READ (SWEEP, ' ( 

READ (SWEEP, ' ( 

CH(42)  =  'S' 

CH(43) 

CH(44) 

CH(45) 

CH(46) 

CH(47) 

N=12 

M=47 

DO  30 


'  )  SF 
CH(12 
CH 
CH 
CH 
CH 
CH 
CH 
CH 
CH 
CH 


) 

(13) 

(14) 

(15) 

(16) 

(17) 

(18) 

(19) 

(20) 

(21) 


(F9.3) ' )  SPF 


Al)') 
IX, Al 
2X,A1 
3X,A1 
4X,A1 
5X,A1 
6X,A1 
7X,A1 
8X,A1 
9X,A1 


CH(26) 

')  CH(27) 

')  CH(28) 

')  CH(29) 

')  CH(30) 

')  CH(31) 

')  CH(32) 

')  CH(33) 

')  CH(34) 

')  CH(35) 


-_   I 


_   I 
—   I 


S 

E' 

S' 

S' 

S" 

S' 


Al) ' )  CH(40) 
IX, Al) ' )  CH(41) 


1=1, N 
WRT(I)  =  0 
DO  40  J=l,< 


c 
c 
c 


C 
C 
C 
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IF( (I-1)*4+J.LE.M)  THEN 

WRT(I)  =WRT(I)  + ICHAR(CH( (I-l)*4+J) )*IP(J) 
ENDIF 


40         CONTINUE 
30      CONTINUE 

CALL  IBWRT(IHP,WRT,M) 

RETURN 

end 


subroutine  error 

common  /ibglob/  ibsta,  iberr,  ibcnt 
write  (*,100)  ibsta, iberr , ibcnt 
100     format  ('  Error ' ,i6, i6,i6) 
return 
end 


SUBROUTINE  TEMPSET(VS) 

common  /ibglob/  ibsta, iberr , ibcnt 

' COMMON  GROUP  1 . ' 

integer*4      cmd(10) ,rd(512) ,wrt(512) ,IP(4) 

CHARACTER* 1  CH(20) 

CHARACTER* 6   VOLTS 

character* 8   bname,bdname 

' COMMON  GROUP  2 . ■ 

LK  =  IBFIND  ('LOCKIN  ') 
IP(4)  -  2**24 
IP(3)  =  2**16 
IP(2)  =  2**8 
IP(1)  =  1 

WRITE ( VOLTS, ' (F6.3) ' )  VS 
CH(1)  =  'X' 
CH(2)  =  '6' 
CH(3)  =  ',' 

READ ( VOLTS, ' (Al) ' )  CH ( 4 ) 
READ(VOLTS, ' (1X,A1) ' )  CH(5) 
READ(VOLTS, ' (2X,A1) ' )  CH(6) 
READ ( VOLTS, ' (3X,A1) ' )  CH(7) 
READ ( VOLTS, ' (4X,A1) ' )  CH(8) 
READ ( VOLTS, ' (5X,A1) ' )  CH(9) 
DO  10  1=1,3 
WRT(I)  =  0 
DO  20  J=l,4 

K=(I-1)*4+J 
IF(K.LE.9)  THEN 

WRT(I)  -  WRT(I)  +  ICHAR(CH(K) )*IP(J) 
ENDIF 
20  CONTINUE 

10       CONTINUE 
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C 
C 
C 


C 
C 
C 


c 
c 
c 


c 
c 
c 


CALL  IBWRT(LK,WRT,9) 

RETURN 

end 


SUBROUTINE  GAINSET(IG) 

common  /ibglob/  ibsta,iberr,ibcnt 

' COMMON  GROUP  1 . ' 

integer*4      cmd(10) ,rd(512) ,wrt(512) ,IP(4) 

CHARACTER* 1  CH(20) 

CHARACTER* 2   GAIN 

character* 8   bname,bdname 

' COMMON  GROUP  2 . ■ 

LK  =  IBFIND  ('LOCKIN  ') 

IP(4)  =  2**24 

IP(3)  =  2**16 

IP(2)  =  2**8 

IP(1)  =  1 

WRITE (GAIN, ' (12) ' )  IG 

CH(1)  =  'G' 

READ(GAIN, ' (Al) ' )  CH(2) 

READ(GAIN, ' (IX, Al) ■ )  CH(3) 

WRT(1)=ICHAR(CH(1) ) +ICHAR(CH( 2 ) )* 
CIP(2)+ICHAR(CH(3) )*IP(3) 

CALL  IBWRT(LK,WRT,3) 

RETURN 

end 


SUBROUTINE  PRES(P) 

character* 14  a 
300     open ( 10 , f ile= ' coml ■ , status= ' old ' ) 

read(10,100)  a 

READ(a, • (3X,F8.3) ',ERR=200)  P 
100     format (al4) 

RETURN 
200      CLOSE (10) 

GOTO  300 

end 


SUBROUTINE  MAX (F,W,FO,M) 
FO=W/FLOAT(M) *FO+F-W/2.0 
RETURN 
END 


SUBROUTINE  KEITH (MV) 
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common  /ibglob/  ibsta,iberr,ibcnt 

' COMMON  GROUP  1 . ' 

integer*4      cmd(10) ,rd(512) ,wrt(512) 

character* 8   bname , bdname , TEMP 

■  COMMON  GROUP  2 . ' 

kl  =  IBFIND  ( 'kill    ' ) 
i24  =  2**24 
il6  =  2**16 
i8  =  2**8 
1006    call  ibrd  (KL,RD,6) 
il  =  rd(l)/i24 
i2  =  (rd(l)-il*i24)/il6 
i3  =  (rd(l)-il*i24-i2*il6)/i8 
i4  =  rd(l)-il*i24-i2*il6-i3*i8 
i5  =  rd(2)/i24 
i6  =  (rd(2)-i5*24)/il6 
i7  =  (rd(2)-i5*i24-i6*il6)/i8 
i88  =  rd(2)-i5*i24-i6*il6-i7*i8 

WRITE (TEMP, ' (8A1) ' )  CHAR(I4) ,CHAR(I3) ,CHAR(I2) , 
CCHAR(Il) ,CHAR(I88) ,CHAR(I7) ,CHAR(I6) ,CHAR(I5) 
READ (TEMP, ' (16) ' ,ERR=1006)  MV 
RETURN 
end 
c 
c 
c 

SUBROUTINE  PLOT ( N , DX , DY , XMAX , XMIN , YMAX , YMIN , NPLOT ) 
REAL*4  DX(1) ,DY(1) 
REAL*4  XMAX , XMIN , YMAX , YMIN 
INTEGER*  2  SYMBOL , PEN , M , DEVTYPE , ADDR 
INTEGER*2  XAXIS,YAXIS 
LOGICAL*2  LIN 

CHARACTER* 40  XLBL,YLBL, TITLE 
C****** **************************************** ********* 

C  GRAPHING  PERAMETERS 

C  eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee 

M=16 
C  eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee 
C  1  IS  FOR  SCREEN  2  IS  FOR  PLOTTER 

DEVTYPE=1 

INTERFACE=#C400 

ADDR=5 
C**************************************************** 

XAXIS=1 

YAXIS=0 

YLBL= '  AMP | ' 

XLBL= '  INDEX  VARIABLE  j  ' 

TITLE= '  AMP  VS  INDEX | ' 

LIN='0' 

SYMBOL=0 

PEN=10 

LEFT=0 . 0 
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BOTTOM=0 . 0 

RIGHT=1.0 

TOP=0 . 8 

IF(NPLOT.EQ.l)  THEN 
CALL  SCREENMODE  (M) 

CALL  PLOTTER  (DEVTYPE, INTERFACE ,ADDR) 
CALL  SCALE  ( XMIN , YMIN , XMAX , YMAX ) 
CALL  WINDOW  (LEFT, BOTTOM, RIGHT, TOP) 

CALL  GRAPHXY ( DX , DY , N , XAXIS , XLBL , YAXIS , 
CYLBL , TITLE , LIN , SYMBOL , PEN ) 

ELSE 

CALL  GRAPHXY (DX,DY,N, XAXIS, XLBL, YAXIS, 
CYLBL , TITLE , LIN , SYMBOL , PEN ) 

ENDIF 

RETURN 

END 


PROGRAM  MAX 3 


C 
C 
C 


C 
C 


100 


INTERFACE 
INTEGER* 4 
END 


TO  SUBROUTINE  TLOOK[C](T) 
T  [NEAR, REFERENCE] 


COMMON  11,18,116,124 

INTEGER*4  11,18,116,124 

INTEGER* 4  NDELY,T,T0 ,TB,TL 

INTEGER  AMP, COUNT, FREQ,FSHIFT 

REAL* 8  ASUM,PHSUM 

REAL* 8  AGSUM 

INTEGER*2  AM( 120 ) ,PH ( 120 ) ,PR( 120 ) 

INTEGER*2  FM( 120 ) ,TM( 120 ), LARGE 

INTEGER* 2  AG (120) 

CHARACTER*20  FILERAW, FILEFORM 

N=120 

NDELY=1000 

11=1 

18=2**8 

116=2**16 

124=2**24 


THE  DATA  OUTPUT  FILE  NAME?' 
THE  FORMAT  OUTPUT  FILE  NAME? 


Write (*,*)  '  WHAT  IS 

READ (*, 100)  FILERAW 

Write (*,*)  '  WHAT  IS 

READ (*, 100)  FILEFORM 

FORMAT (A) 

open ( 3 , f ile= ' block ' , status= ' new ' ) 

OPEN ( 5 , FILE=FILEFORM, STATUS=  *  NEW ' ) 

OPEN (10,FILE=FILERAW,FORM=* UNFORMATTED'  , STATUS =' NEW '  ) 

write (*,*)  '  WHAT  IS  THE  RESONANCE  DAC  SETTING?' 

READ(*,*)  ICENT 
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SWEEP? 


c 
c 

c 


200 


300 


400 


C 

c 
c 
c 
c 
c 


700 


500 


1000 


ISTART=ICENT-N/2 

CALL  VOLTS0(IFLIP(ISTART) ) 

Write(*,*)  *  WHAT  IS  THE  INITIAL  TEMPERATURE  VOLTAGE? 

READ(*,*)  TBB 

ITB=INT(TBB) 

CALL  VOLTS 1 ( IFLIP ( ITB ) ) 

Write (*,*)  '  WHAT  IS  THE  FINAL  TEMPERATURE  VOLTAGE? 

READ(*,*)  TF 


Write(*,*) 


HOW  MANY  HOURS  FOR  THIS  TEMPERATURE 


READ(*,*)  TIM 
TIM=TIM*3600000.0 
TMM=(TF-TBB)/TIM 
CALL  TLOOK(T) 
T0=T 
WRITE (*,*)  '  WHAT  IS  THE  SENSITIVITY  ID  NUMBER? 

READ(*,*)  IG 

OPEN ( 15 , FILE= ' SENSE ' , STATUS= ' NEW ' ) 
CALL  TEMPGET(TEMP) 
TSHIFT=TEMP-22.0 
ITSHIFT=TSHIFT 
WRITE (5, 200)  TSHIFT 

FORMAT ( IX ,  ' TEMPERATURE  SHIFT= ' , F8 . 3 ) 
CALL  PRESGET(P) 
PSHIFT=P-3000.0 
WRITE (5, 300)  PSHIFT 

FORMAT( IX, 'PRESSURE  SHIFT= ' , F10 . 3 ) 
CALL  FREQGET ( FREQ ) 
FSHIFT=FREQ-10000 
WRITE (5, 400)  FSHIFT 
FORMAT (IX, 'FREQUENCY  SHIFT=',I6) 
CALL  PHASEGET(AMPL) 

CALL  LOCKIN(AMPL,ANGL) 

ASHIFT=ANGL-80.0 

WRITE (5, 7 00)  ASHIFT 

FORMAT ( IX, 'ANGLE  SHIFT= ' ,F8 . 3 ) 

WRITE ( * , * )  ' AMPL= ' , AMPL 

WRITE ( * , * )  ' ANGL= ' , ANGL 
WRITE (5, 500)  N 
FORMAT ( IX , '  LOOP  RANGE= ' , 14 ) 
CLOSE ( 5 ) 
KBLOCK=0 
KBLOCK=KBLOCK+ 1 
LARGE=-20000 
COUNT=0 

ISTART=ICENT-N/2 
ISTOP=ISTART+N 
CALL  VOLTS0 ( IFLIP ( ISTART ) ) 
CALL  SECS(T,T0) 
TEMP=FLOAT ( T ) *TMM+TBB 
ITEMP=INT(TEMP) 
CALL  VOLTS 1( IFLIP ( ITEMP ) ) 
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C . . . FORWARD  SWEEP 

DO  20  I=ISTART,ISTOP 

CALL  VOLTS0 ( IFLIP ( I ) ) 
CALL  SECS(T,T0) 
TB=T 

ASUM=0 . 0 
PHSUM=0 . 0 
C  AGSUM=0 . 0 

KOUNT=0 
3000  CALL  SECS(T,T0) 

CALL  Kel77b(AMP) 
C  CALL  PHASEGET(AMPL,ANGL) 

CALL  PHASEGET(AMPL) 
ASUM=ASUM+FLOAT ( AMP ) 
PHSUM=PHSUM+AMPL 
C  AGSUM=AGSUM+ANGL 

C  WRITE ( * , * )  ' PHSUM= ' , PHSUM 

C  WRITE ( * , * )  ' ASUM= ' , ASUM 

C  WRITE ( * , * )  ' AGSUM= ' , AGSUM 

KOUNT=KOUNT+l 
TL=T-TB 
IF(TL.LT.NDELY)  GOTO  3000 
CALL  TEMPGET(TEMP) 
CALL  FREQGET ( FREQ ) 
CALL  PRESGET(P) 
COUNT=COUNT+l 

AM ( COUNT )=INT2(ASUM/FLOAT(KOUNT) ) 
PH ( COUNT ) =INT2 ( ( PHSUM/FLOAT ( KOUNT ) ) * 1 . 0E6 ) 
C  AG ( COUNT ) =INT2 ( ( ( AGSUM/FLOAT ( KOUNT ) ) 

C      C-ASHIFT)*100.0) 

C  WRITE  (*,*)  ,AM=' , AM (COUNT) 

C  WRITE  (*,*)  'PH=' ,PH( COUNT) 

C  WRITE  (*,*)  'AG=' , AG (COUNT) 

FM ( COUNT ) =INT2 ( FREQ-FSHIFT ) 
TM( COUNT )=INT2( (TEMP-TSHIFT) *1000 . 0 ) 
PR ( COUNT )=INT2( (P-PSHIFT) *10 . 0 ) 
LARGE=MAX0 ( LARGE , AM ( COUNT ) ) 
IF  ( LARGE. EQ. AM ( COUNT ) )  ICENT=I 
C  CALL  ADJUST (IG) 

C  ISEN( COUNT )=INT1(IG) 

20       CONTINUE 
C         WRITE(IO)  AM,PH,PR,FM,TM,AG 

WRITE (10)  AM,PH,PR,FM,TM 
C         CALL  SENSET(IG) 

WRITE (*, 30)  AM(120) ,PH(120) ,PR(120) ,FM(120) , 
CTM(120)+ITSHIFT*1000 
C         WRITE (15,*)  IG 
C      C,AG(120) 
30      FORMAT (IX, 51 10) 
CALL  SECS(T,T0) 
IF(T.GT.INT(TIM) )  GOTO  2000 
Write ( * , * )  ' BLOCK  NUMBER ■ , KBLOCK 
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GOTO  1000 
2000    Write (3, 600)  KBLOCK 
600     FORMAT(lX,'  NUMBER  OF  BLOCKS  OF  DATA=',I5) 

CLOSE (10) 

CLOSE ( 3 ) 
s        CLOSE (15) 

END 


C 
C 
C 


C 
C 
C 


C 
C 
C 


C 
C 
C 


FUNCTION  IFLIP( INPUT) 

IFLIP=-INPUT+4095 

RETURN 

END 


SUBROUTINE  VOLTS 0 ( ITEMP ) 

INTEGER*  2  ADAPT , DEVICE , CTRL , STAT , VDACH 1 

VDACH1  =  INT2( ITEMP) 

ADAPT=0 

DEVICE  =  9 

CTRL  =  0 

STAT  =  0 

CALL  AOUS  ( ADAPT, DEVICE, 0, CTRL, VDACH 1, STAT) 

RETURN 

END 


SUBROUTINE  VOLTS 1 ( ITEMP ) 

INTEGER* 2  ADAPT, DEVICE , CTRL, STAT, VDACH 1 

VDACH 1  =  INT2( ITEMP) 

ADAPT=0 

DEVICE  =  9 

CTRL  =  0 

STAT  =  0 

CALL  AOUS  ( ADAPT, DEVICE, 1, CTRL, VDACH 1, STAT) 

RETURN 

END 


SUBROUTINE  PHASEGET(AMPL) 

common  /ibglob/  ibsta,iberr,ibcnt 

INTEGER*4  11,188,1166,124 

DIMENSION  1(20) 

■  COMMON  GROUP  1 .  ' 

integer  CR 

integer  cmd(10) ,rd(512) ,wrt(512) ,IP(4) 

CHARACTER* 1  CH(20) 

character* 8   bname,bdname 
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C 
C 

c 

C  1004 
C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

cc 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 
c 
c 
c 
c 
c 
c 
c 
c 


CHARACTER* 20  TEMPO, TEMP 
' COMMON  GROUP  2 . ■ 
data  CR/13/ 

LK  =  IE-FIND  ('LOCKIN  ') 
IP(4)  =  2**24 
IP(3)  =  2**16 
IP(2)  -  2**8 
IP(1)  =  1 
124  =  2**24 
1166  =  2**16 
188   =  2**8 
11=1 
CH(1)  =  'P' 
WRT(1)=ICHAR(CH(1) ) 
CALL  IBWRT(LK,WRT,1) 
IBRD  (LK,RD,12) 
=  RD(1)/I24 

=  (RD(1)-I(4)*I24)/I16 
=  (RD(1)-I(4)*I24-I(3) 
=  RD(1)-I(4)*I24-I(3) 
=  RD(2)/I24 

=  (RD(2)-I(8)*I24)/I16 
=  (RD(2)-I(8)*I24-I(7) 
RD(2)-I(8)*I24-I(7) 
RD(3)/I24 

(RD(3)-I(12)*I24)/I 
(RD(3)-I(12)*I24-I( 
RD(3)-I(12)*I24-I(1 


10 
30 


CALL 
1(4) 
1(3) 
1(2) 

1(1) 
1(8) 

1(7) 

1(6) 

1(5) 

1(12)  = 

1(11)  = 

1(10)  = 

1(9)  = 

DO  10  J=l,12 

WRITE  (*,*)  I (J) 

IF  (I(J) .EQ.CR) 
JJ=J-1 
GOTO  30 

END  IF 
CONTINUE 
IF  (JJ.EQ.7)  THEN 

WRITE (TEMP, ' (7A1) 


*I166)/I88 
*I166-I(2)*I88 


,  CHAR(I( 
THEN 


*I166)/I88 
*I166-I(6)*I88 

166 

11)*I166)/I88 

1)*I166-I(10)*I88 

J)) 


CCHAR(I(3)),  CCHAR(I 

CCHAR(I(6) ) ,CHAR(I(7 

READ ( TEMP , ' ( F7 . 2 

ELSE  IF  (JJ.EQ.6 


WRITE(TEMP, ' (6A1) ' )  CHAR(I(1) ) ,CHAR(I(2) ) , 

)),CHAR(I(6)) 
)  ANGL 


CCHAR(I(3) ) ,CHAR(I(4 
READ ( TEMP , ' ( F6 . 2 
ELSE  IF  (JJ.EQ.5 


CCHAR(I(3) ) ,CHAR(I(4 
READ ( TEMP , ' ( F5 . 2 
ELSE  IF  (JJ.EQ.4 
WRITE (TEMP, ' (4A1) 

CCHAR (1(3)), CHAR (1(4 
READ(TEMP, ' (F4.2 


)  CHAR(I(1) ) ,CHAR(I(2) ) , 
(5)), 


4) ),CHAR(I 

) 

•,ERR=1004)  ANGL 

THEN 


),CHAR(I(5 
,ERR=1004 
THEN 


WRITE(TEMP, ' (5A1) ' )  CHAR(I(1) ) ,CHAR(I(2) ) , 

)) 

)  ANGL 


) , CHAR (1(5 
,ERR=1004 
THEN 
')  CHAR(I 


(1)),CHAR(I(2)), 
,ERR=1004)  ANGL 
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C  END  IF 

CH(1)  =  'Q' 
WRT(1)=ICHAR(CH(1) ) 
CALL  IBWRT(LK,WRT,1) 
1005     CALL  IBRD  (LK,RD,12) 
1(4)  =  RD(1)/I24 
1(3)  =  (RD(1)-I(4)*I24)/I166 
1(2)  =  (RD(1)-I(4)*I24-I(3)*I166)/I88 
1(1)  =  RD(1)-I(4)*I24-I(3)*I166-I(2)*I88 
1(8)  =  RD(2)/I24 
1(7)  =  (RD(2)-I(8)*I24)/I166 
1(6)  =  (RD(2)-I(8)*I24-I(7)*I166)/I88 
1(5)  =  RD(2)-I(8)*I24-I(7)*I166-I(6)*I88 
1(12)  =  RD(3)/I24 
1(11)  =  (RD(3)-I(12)*I24)/I166 
1(10)  -  (RD(3)-I(12)*I24-I(11)*I166)/I88 
1(9)  =  RD(3)-I(12)*I24-I(11)*I166-I(10)*I88 
DO  20  K=l,12 
IF  (I(K).EQ.CR)  THEN 
KK=K-1 
GOTO  40 
END  IF 
20      CONTINUE 
40      IF  (KK.EQ.10)  THEN 

WRITE (TEMPO, ' (10A1) ' )  CHAR(I(1) ),CHAR(I(2) ), 
CCHAR(I(3) ) ,CCHAR(I(4) ),CHAR(I(5) ) , 
CCHAR(I(6) ) ,CHAR(I(7) ),CHAR(I(8) ) , 
CCHAR(I(9) ) ,CHAR(I(10)) 

READ (TEMPO, ' (F10.2) ',ERR=1005)  AMPL 
C  IF  (1(10). EQ. 51)  AMPL=AMPL/ 10. 0**3.0 

C  IF  (1(10) .EQ. 54)  AMPL=AMPL/ 10. 0**6.0 

C  IF  (1(10). EQ. 57)  AMPL=AMPL/ 10. 0**9.0 

ELSE  IF  (KK.EQ.9)  THEN 
WRITE(TEMPO, ' (9A1) ' )  CHAR(I(1) ) ,CHAR(I(2) ) , 
CCHAR(I(3) ) ,  CCHAR(I(4) ),CHAR(I(5) ) , 
CCHAR(I(6) ) ,CHAR(I(7) ) ,CHAR(I(8) ),CHAR(I(9) ) 
READ(TEMPO, ' (F9.2) ' ,ERR=1005)  AMPL 
C  IF  (I(9).EQ.51)  AMPL=AMPL/ 10. 0**3.0 

C  IF  (I(9).EQ.54)  AMPL=AMPL/ 10. 0**6.0 

C  IF  (I(9).EQ.57)  AMPL=AMPL/ 10. 0**9.0 

ELSE  IF  (KK.EQ.8)  THEN 
WRITE ( TEMPO, • ( 8A1 ) ■ )  CHAR (1(1)), CHAR (1(2)), 
CCHAR(I(3) ) ,  CCHAR(I(4) ) ,CHAR(I(5) ) , 
CCHAR(I(6) ) ,CHAR(I(7) ) ,CHAR(I(8) ) 

READ(TEMPO, ' (F8.2) ' ,ERR=1005)  AMPL 
C  IF  (I(8).EQ.51)  AMPL=AMPL/10. 0**3.0 

C  IF  (I(8).EQ.54)  AMPL=AMPL/ 10. 0**6.0 

C  IF  (I(8).EQ.57)  AMPL=AMPL/ 10. 0**9.0 

ELSE  IF  (KK.EQ.7)  THEN 
WRITE (TEMPO, ' (7A1) ')  CHAR(I(1) ),CHAR(I(2) ), 
CCHAR(I(3)),  CCHAR(I(4)),CHAR(I(5)),CHAR(I(6)),CHAR(I(7)) 
READ(TEMPO, ' (F7.2) ' ,ERR=1005)  AMPL 
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C  IF  (I(7).EQ.51)  AMPL=AMPL/ 10. 0**3.0 

C  IF  (I(7).EQ.54)  AMPL=AMPL/10. 0**6.0 

C  IF  (I(7).EQ.57)  AMPL=AMPL/ 10. 0**9.0 

END  IF 

IF  (ABS(AMPL) .GT. 20000.0)  GOTO  1005 

RETURN 

end 


C 
C 

c 


SUBROUTINE  Kel77b(MV) 

common  /ibglob/  ibsta, iberr, ibcnt 

COMMON  11,18,116,124 

INTEGER*4  11,18,116,124 

' COMMON  GROUP  1 . ■ 

integer*4      cmd(10) ,rd(512) ,wrt(512) 

character* 8   bname , bdname , TEMP 

' COMMON  GROUP  2 . ' 

K177BAD=IBFIND( "K177B  ') 
1006    call  ibrd  (K177BAD,RD, 6 ) 

il  =  rd(l)/i24 

i2  =  (rd(l)-il*i24)/il6 

i3  =  (rd(l)-il*i24-i2*il6)/i8 

i4  =  rd(l)-il*i24-i2*il6-i3*i8 

15  =  rd(2)/i24 

i6  =  (rd(2)-i5*i24)/il6 

i7  =  (rd(2)-i5*i24-i6*il6)/i8 

i88  =  rd(2)-i5*i24-i6*il6-i7*i8 

WRITE (TEMP, ' (8A1) ■ )  CHAR(I4) ,CHAR(I3) ,CHAR(I2) 
CCHAR(Il) ,CHAR(I88) ,CHAR(I7) ,CHAR(I6) ,CHAR(I5) 

READ (TEMP, * (16) ' ,ERR=1006)  MV 

IF  (ABS(MV) .GT. 20000)  GOTO  1006 
RETURN 
end 


C 
C 
C 


SUBROUTINE  FREQGET ( FREQ ) 

common  /ibglob/  ibsta, iberr, ibcnt 

COMMON  11,18,116,124 

INTEGER*4  11,18,116,124 

INTEGER  FREQ 

' COMMON  GROUP  1 . ' 

integer  SENAD 

integer*4      cmd(10) ,rd(512) ,wrt(512) 

character* 8   bname, bdname 

CHARACTER* 7  TEMP 

' COMMON  GROUP  2 . ' 

i3  =  ichar( '3' ) 

iO  =  ichar( '0' ) 

SENAD=IBFIND( ' SENCORE  ') 
1006    call  ibrd  (SENAD, RD, 17 ) 
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C 
C 
C 


100 


IC(1)  =  rd(l)/i24 

IC(2)  =  (rd(l)-IC(l)*i24)/il6 

IC(3)  =  (rd(l)-IC(l)*i24-IC(2)*il6)/i8 

IC(4)  =  rd(l)-IC(l)*i24-IC(2)*il6-IC(3)*i8 

IC(5)  =  rd(2)/i24 

IC(6)  =  (rd(2)-IC(5)*i24)/il6 

IC(7)  =  (rd(2)-IC(5)*i24-IC(6)*il6)/i8 

IC(8)  =  rd(2)-IC(5)*i24-IC(6)*il6-IC(7)*i8 

IC(9)  =  rd(3)/i24 

IC(10)  =  (rd(3)-IC(9)*i24)/il6 

IC(ll)  =  (rd(3)-IC(9)*i24-IC(10)*il6)/i8 

IC(12)  =  rd(3)-IC(9)*i24-IC(10)*il6-IC(ll)*i8 

IC(13)  =  rd(4)/i24 

IC(14)  =  (rd(4)-IC(13)*i24)/il6 

IC(15)  =  (rd(4)-IC(13)*i24-IC(14)*il6)/i8 

IC(16)  =  rd(4)-IC(13)*i24-IC(14)*il6-IC(15)*i8 

IF(IC(14) .EQ.i3)  THEN 

WRITE ( TEMP , ' ( 7A1 ) ' )  CHAR ( IC ( 1 ) ) , CHAR ( IC ( 8 ) ) , 

CCHAR(IC(7) ) ,  CHAR(IC(5) ) ,CHAR(IC(12) ), 

CCHAR(IC(11) ) ,CHAR(IC(10) ) 

READ ( TEMP, ' (17) ',ERR=1006)  FREQ 

ENDIF 

IF(IC(14) .EQ.iO)  THEN 

WRITE(TEMP, ' (7A1) ■ )  CHAR(IC(1) ) ,CHAR(IC(8) ) , 

CCHAR(IC(7) ) ,CHAR(IC(6) ) ,CHAR(IC(5) ), 

CCHAR(IC(12) ) ,CHAR(IC(10) ) 

READ (TEMP, ■ (17) ' ,ERR=1006)  FREQ 
ENDIF 
RETURN 
end 


SUBROUTINE  SECS(T,T0) 

INTEGER* 4  T, TO, DAY 

DAY=3600*24*1000 

CALL  TLOOK(T) 

IF(T.GT.DAY.OR.T.LT.O)  GOTO  100 

T=T-T0 

IF(T.LT.O)  THEN 

T0=T0-DAY 

T=T+DAY 
ENDIF 
RETURN 
END 


C 
C 
C 


SUBROUTINE  TEMPGET(T) 

common  /ibglob/  ibsta, iberr,ibcnt 

COMMON  11,18,116,124 

INTEGER*4  11,18,116,124 
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' COMMON  GROUP  1 . ' 

integer*4      cmd(10) ,rd(512) ,wrt(512) 

character* 8   bname , bdname , TEMP 

' COMMON  GROUP  2 . ' 

K195AAD=IBFIND( 'K195A  ') 
1004     CALL  IBRD  (K195AAD,RD, 12 ) 

15  =  RD(2)/I24 

16  =  (RD(2)-I5*I24)/I16 

17  =  (RD(2)-I5*I24-I6*I16)/I8 
188  =  RD(2)-I5*I24-I6*I16-I7*I8 
19  =  RD(3)/I24 

110  =  (RD(3)-I9*I24)/I16 

111  =  (RD(3)-I9*I24-I10*I16)/I8 

112  =  RD(3)-I9*I24-I10*I16-I11*I8 
WRITE(TEMP, ' (8A1) ' ) char ( 188 ), char ( 17 ) , 

Cchar(l6) ,char(I5) ,  CCHAR(I12) ,CHAR(I11) , 
CCHAR(I10),CHAR(I9) 

READ(TEMP, ' (F8.3) ' ,ERR=1004)  R 

R0=99.98 

ALPHA=0. 0039076 

DELTA=1.5205 

ALDEL=ALPHA* DELTA 

RC=R/R0-1.0 

PARTl=ALPHA+ALDEL/100 . 0 

PART2=(ALDEL/100.0+ALPHA)**2 

PART3=4 . 0*RC*ALDEL/ 10000 . 0 

PART4=2 . 0*ALDEL/10000 . 0 

TEMP2=PART1-SQRT ( PART2-PART3 ) 

T=TEMP2/PART4 

RETURN 

end 


c 
c 


subroutine  error 

common  /ibglob/  ibsta,  iberr,  ibcnt 
write  (*,100)  ibsta, iberr, ibcnt 
100    format  ('  Error ■ ,i6,i6, i6 ) 
return 
end 


C 
C 
C 


SUBROUTINE  PRESGET(P) 

common  /ibglob/  ibsta, iberr, ibcnt 

COMMON  11,18,116,124 

INTEGER*4  11,18,116,124 

' COMMON  GROUP  1 . ' 

integer  BMAD 

integer*4      cmd( 10) ,rd(512) ,wrt (512) 

character* 8  bname, bdname 

character*7   Pres 
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' COMMON  GROUP  2 . ' 

BMAD=IBFIND( 'BM600  ') 
1006     call  ibrd  (BMAD,RD,10) 

il  =  rd(l)/i24 

i2  =  (rd(l)-il*i24)/il6 

i3  =  (rd(l)-il*i24-i2*il6)/i8 

i4  =  rd(l)-il*i24-i2*il6-i3*i8 

i5  =  rd(2)/i24 

i6  =  (rd(2)-i5*i24)/il6 

i7  =  (rd(2)-i5*i24-i6*il6)/i8 

i88  =  rd(2)-i5*i24-i6*il6-i7*i8 

WRITE(Pres, ■ ( 7A1 ) ■ )  CHAR(I4) ,CHAR(I3) ,CHAR(I2) 
CCHAR(Il) ,CHAR(I88) ,CHAR(I7) ,CHAR(I6) 

READ(PRES, ' (F7.1) ' ,ERR=1006)  P 
RETURN 
end 


C 

C 


SUBROUTINE  SENSET(IG) 

common  /ibglob/  ibsta,iberr,ibcnt 

INTEGER*4  II , 188 , 1166 , 124 

DIMENSION  1(20) 

' COMMON  GROUP  1 . ' 

integer  CR 

integer  cmd(10) ,rd(512) ,wrt(512) ,IP(4) 

CHARACTER* 1  CH(20) 

CHARACTER* 2  GAIN 

character* 8   bname,bdname 

CHARACTER* 20  TEMP 

' COMMON  GROUP  2 . ' 

data  CR/13/ 

IP(4)  =  2**24 

IP(3)  =  2**16 

IP(2)  =  2**8 

IP(1)  =  1 

124  =  2**24 

1166  =  2**16 

188   =  2**8 

11=1 

J=0 

K=l 

L=0 
20      LK  =  IBFIND  ('LOCKIN  ') 

CH(1)  =  'Y' 

CH(2)  =  '4' 

WRT(1)=ICHAR(CH(1) ) +ICHAR(CH( 2 ) )*IP(2) 

CALL  IBWRT(LK,WRT,2) 
1004     CALL  IBRD  (LK,RD,12) 

1(4)  =  RD(1)/I24 

1(3)  =  (RD(1)-I(4)*I24)/I166 

1(2)  =  (RD(1)-I(4)*I24-I(3)*I166)/I88 

1(1)  =  RD(1)-I(4)*I24-I(3)*I166-I(2)*I88 
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WRITE (TEMP, ■ (Al) ' )  CHAR(I(1)) 

READ (TEMP, ' (II) ' ,ERR=1004)  ISTATUS 

L=L+1 

IF  (L.LT.3)  GOTO  20 

IF  ( J.  EQ.  LAND.  ISTATUS.  EQ.O)  GOTO  30 

IF  (J. EQ.l. AND. K. EQ.O. AND. ISTATUS. EQ.O)  GOTO  30 

IF  (ISTATUS.EQ.l)  THEN 

IG=IG+2 

IF  (IG.GT.24)  IG=24 

IF  (IG.LT.l)  IG=1 

J=l 

GOTO  10 
ELSE  IF  (ISTATUS. EQ.O)  THEN 

IG=IG-1 

K=0 

IF  (IG.GT.24)  IG=24 

IF  (IG.LT.l)  IG=1 

GOTO  10 
END  IF 
10       WRITE (GAIN, ' (12) ' )  IG 
CH(1)  =  -G' 

READ(GAIN, ' (Al) ' )  CH(2) 
READ(GAIN, ' (IX, Al) ■ )  CH(3) 

WRT(1)=ICHAR(CH(1) ) +ICHAR(CH ( 2 ) )*IP(2)+ 
CICHAR(CH(3) )*IP(3) 

CALL  IBWRT(LK,WRT,3) 
GOTO  20 
30       RETURN 
end 


PROGRAM  VIKAN2 

$large 

INTERFACE  TO  SUBROUTINE  TIME  (N,STR) 

CHARACTER* 10  STR  [NEAR, REFERENCE] 

INTEGER* 2  N  [VALUE] 

END 

REAL*4  DX(3000) ,DY(3000) ,XMAX,XMIN, YMAX, YMIN 

CHARACTER* 1  KOL 

CHARACTER* 2  SWEEP 

CHARACTER* 10  DNAME 

CHARACTER* 7  DTEMP 

CHARACTER* 10  TSTR 

INTEGER*2  HR,MIN,SEC 

INTEGER  READS, T,TC, TO, TOO, POINTS 

SWEEP='15' 

WIDTH  =30.0 

READS  =  3 

KOUNT  =  0 

STP  =  -0.00125 

CALL  KEITHA(TZERO) 
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OPEN(60,File='c:DATA.BUF' , STATUS= ' OLD • ) 

READ (60,*)  VOLT,IG 

CALL  GAINSET(IG) 

CLOSE (60) 

VOLTE  =  (10.0-TZERO)*0.01  +  VOLT 

write  (*,*)  volte 

CLOSE (10) 

DTEMP='c:DATA. ' 

OPEN (70, File= ■ c : NAME . NUM ' , STATUS= ■ OLD ■ ) 

READ(70,301)  DNAME 
301      FORMAT ( IX, A10) 

CLOSE  (70) 

READ ( DNAME, ' (7X, 13) * )  NUM 

NUM=NUM+1 

OPEN (70, File=  *  c : NAME . NUM ' , STATUS= ' NEW • ) 

WRITE (70, 302)  DTEMP,NUM 

WRITE (*, 302)  DTEMP,NUM 
302      FORMAT ( IX, A7, 13) 
CLOSE  (70) 

OPEN (20, FILE=DNAME , STATUS= ■ NEW • ) 

OPEN(30,File='c:Peak.buf ' , STATUS= ' OLD • ) 

READ(30,*)  F,TI 

TI=TI+273.15 

CALL  KEITHA(TF) 

TF=TF+273.15 

F=SQRT(TF/TI)*F 

CLOSE (30) 

IXT=3 
c  IXT=1 

NA  =  21/IXT 

IDELY=2*NA 

Q*  ********************************************************  * 

c  M=1000/NA/IXT-IDELY/NA 

C  USED  ONLY  FOR  CENTERING 

C********************************************************** 

M=29 

C  USED  AFTER  CENTERING 
C********************************************************** 

XMAX=FLOAT(M) 
1001     CALL  TIME(10,TSTR) 

READ(TSTR, ' (12) ' ,ERR=1001)  HR 
READ(TSTR, ' (3X,I2) ' ,ERR=1001)  MIN 
READ(TSTR, ' (6X,I2) ' ,ERR=1001)  SEC 
TOO  =  HR*3600+MIN*60+SEC 

100  CALL  TEMPSET(VOLT) 

101  DO  30  1=1, READS 
500        SF=F-WIDTH/2.0 

SPF=F+WIDTH/2.0 
INDEX=1 
POINTS=M 
ISTART=1 
200        CALL  HP (SF,SPF, SWEEP) 


185 


AMAX=0 . 0 
XMIN=200000.0 
DO  42  ID=1,IDELY 
CALL  KEITH (MV) 
42         CONTINUE 

DO  40  J=ISTART, POINTS, INDEX 
AVG=0 . 0 
DO  41  JA=1,NA 
300  CALL  KEITH (MV) 

KOUNT=KOUNT+l 
IF(KOUNT.EQ.2*NA)  THEN 
IG=IG+1 

CALL  GAINSET(IG) 
ENDIF 

IF(KOUNT.GT.3*NA)  THEN 
KOUNT=0 
GOTO  500 
ENDIF 

IF(MV. GT.20000.OR.MV.LT. 100)  GOTO  300 
KOUNT=0 

AVG=AVG+FLOAT ( MV ) 
41         CONTINUE 

AMP=AVG / FLOAT ( NA ) 
DX(J)=FLOAT(J) 
DY ( J ) =AMP 
IF(AMP.GT.AMAX)  THEN 
AMAX=AMP 
JMAX=J 
ENDIF 

IF(AMP.LT.XMIN)  XMIN=AMP 
40         CONTINUE 

IF ( INDEX. EQ.l)  THEN 

NPLOT=l 
ELSE 

NPLOT=0 
ENDIF 
YMAX=AMAX 
C********************* ************  ********************** 

c  CALL  PLOT(M,DX,DY,XMAX,XMIN,YMAX,YMIN,NPLOT) 

C  ONLY  USED  FOR  CENTERING 
C******************************************************* 

IF (INDEX. EQ.l)  THEN 
INDEX=-INDEX 
TEMP=SF 
SF=SPF 
SPF=TEMP 
ISTART=POINTS 
POINTS=l 
JFMAX=JMAX 
AVGA=AMAX 
CALL  KEITHA(TI) 
IF(TI.LT.IO.O)  GOTO  169 
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CALL  PRES(VI) 
GOTO  200 
ELSE 

CALL  KEITHA(TF) 
CALL  PRES(VF) 
JRMAX=JMAX 

FO=FLOAT ( JFMAX+ JRMAX ) /2 . 0 
CALL  MAX(F, WIDTH, FO,M) 
F=FO 

AMP= ( AMAX+AVGA ) / 2 . 0 

IF ( JFMAX . EQ . 1 . OR . JFMAX . EQ . M . OR . JRMAX . EQ . 1 . OR . 
C JRMAX. EQ.M)  THEN 

OPEN(30,File='c:Peak.buf ' , STATUS= ■ OLD ' ) 
READ(30,*)  F,TI 
CLOSE (30) 
TI=TI+273.15 
TF=TF+273.15 
F=SQRT(TF/TI)*F 
GOTO  500 
ENDIF 
1003  CALL  TIME(10,TSTR) 

READ(TSTR, ' (12) ' ,ERR=1003)  HR 
READ(TSTR, ' (3X,I2) ' ,ERR=1003)  MIN 
READ(TSTR, ' (6X,I2) ' ,ERR=1003)  SEC 
T  =  HR*3600+60*MIN+SEC 
T=T-T00 

IF(T.LT.O.O)  THEN 
T00=T00-24*3600 
T  =  HR*3600+60*MIN+SEC 
T=T-T00 
ENDIF 
IF(AMP.LT.7000)  THEN 
IG=IG-1 

IF(IG.GE.7)  CALL  GAINSET(IG) 
ENDIF 

IF(AMP.GT. 19000)  THEN 
IG=IG+1 

IF(IG.LE.24)  CALL  GAINSET(IG) 
ENDIF 

RATIO=F/SQRT( (TI+TF+273 . 15*2 . 0 ) /2 . 0 ) 
WRITE (20, *)  TI,TF,VI,VF,AMP,F,T,IG 
C*********************************************************** 

c  DO  18  IF=1,19 

c  WRITE ( * , * ) 

c  18          CONTINUE 

c  WRITE(*,390)  SPF,SF 

c  390         FORMAT ( IX, 'FREQUENCY=  ' ,F10 . 3 , 48X,F10 . 3 ) 

c  WRITE (*,*) 

c  IF  (JRMAX  .GT.  JFMAX  )  M=M-1 

c  IF  (JRMAX  .LT.  JFMAX  )  M=M+1 

c  DF=0.0 

C  USED  ONLY  FOR  CENTERING 


187 

Q*********** ******  ******  ************************************ 

WRITE (*,400)TI,TF,VI,VF,M,F,T,JRMAX,JFMAX, RATIO 
400       FORMAT(1X,2F8.3,2F8.1,I6,F10.3,I8,2I4,F9.3) 

Q** ***************************************************  ****** 

DT=TF-TI 
DF=FLOAT ( JRMAX- JFMAX ) /FLOAT ( M ) *WIDTH 
C   USED  ONLY  AFTER  CENTERING 

Q*********************************************************** 

F=F+DF*1.5 
ENDIF 
30    CONTINUE 

IF ( ABS ( DT ) . GT . 0 . 0  2 )  READS=READS+ 1 
IF ( ABS ( DT ) . GT . 0 . 0  3 )  READS=READS+ 1 
IF(ABS(DT) .GT.0.04)  READS=READS+1 
IF(ABS(DT) .GT.0.05)  READS=READS+1 
IF ( ABS ( DT ) . GT . 0 . 0  6 )  READS=READS+ 1 
IF(ABS(DT) .GT.0.07)  READS=READS+1 
IF(ABS(DT) .GT.0.08)  READS=READS+ 1 
C        IF(ABS(DT) .LT. 0.05. AND. READS. GT. 2)  READS=READS-1 
C        IF(ABS(DT) .LT. 0.04. AND. READS. GT. 2)  READS=READS-1 
C        IF(ABS(DT) .LT. 0.03. AND. READS. GT. 2)  READS=READS-1 
IF(ABS(DT) . LT. 0.02. AND. READS. GT. 2)  READS=READS-1 
IF(ABS(DT) . LT. 0.015. AND. READS. GT. 2)  READS=READS-1 
IF(ABS(DT) . LT. 0.01. AND. READS. GT. 2)  READS=READS-1 
IF(ABS(DT) . LT. 0.005. AND. READS. GT. 2)  READS=READS-1 
WRITE (*, 380)  VOLT 
380    FORMAT  (IX,' SET  POINT  VOLTAGE  =  ' ,F9.4) 
c       Call  scan(F,IG) 

OPEN(30,File='c:Peak.buf ■ , STATUS =' NEW ' ) 
WRITE(30,*)  F,TF 
CLOSE (30) 

if  (tf.lt. 10.0)  goto  169 
IF(TF.GT. 100.0)  GOTO  169 
IF ( VOLT. GE. VOLTE)  THEN 
VOLT=VOLT+STP 

OPEN(60,File=,c:VOLT.BUF' , STATUS= ' NEW ' ) 
WRITE (60,*)  VOLT 
CLOSE (60) 
GOTO  100 
ELSE 

IF(VOLT.GT.VOLTE-0.4)  THEN 
VOLT=VOLTE-0.5 

OPEN(60,File='c:VOLT.BUF' , STATUS= ' NEW ' ) 
WRITE (60,*)  VOLT 
CLOSE (60) 
GOTO  100 
ENDIF 
GOTO  101 
ENDIF 
169    OPEN ( 70 ,File='c: STOP. RUN' , STATUS='NEW' ) 
CLOSE (70) 
CLOSE (20) 
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Call  scan(F,IG) 

Call  tempset(O.O) 

STOP 

END 
c 
c 
c 

Subroutine  scan(F,IG) 

CHARACTER* 2  SWEEP 

Character* 3  name 

Character* 8  fname 

Sweep='60' 
c        SWEEP=,15' 

OPEN (80, FILE= ' INDEX . BUF ' , STATUS= ' OLD ' ) 

READ (80,*)  INDEX 

CLOSE (80) 

INDEX=INDEX+1 

OPEN (80, FILE= ' INDEX . BUF ' , STATUS= ■ NEW ' ) 

WRITE (80,*)  INDEX 

CLOSE (80) 

WRITE(NAME, ' (13) ' )  INDEX 

Write  (fname,' (Al,  A3,  A4 ) • )  • F' , NAME, ' .DAT' 

Open  ( 50 , f ile=f name , status= ' new  * ) 

WIDTH  =  100.0 

IF  (IG.LT.24)  IG=IG+1 

CALL  GAINSET(IG) 

WRITE (50,*)  IG,F 
c    the  following  line  is  for  avoiding  too  much  signal 
c    at  the  end  of  the  run. 

f=f+200 

SF  =  F-WIDTH/2.0 

SPF  =  F+WIDTH/2.0 
C        XMAX=SPF 
C        XMIN=SF 
C        YMAX=-1.0E32 
C        YMIN=-YMAX 

IXT=3 
C        IXT=1 

NA=15/IXT 

NP=2790/NA/IXT 

SLOPE=WIDTH/ FLOAT ( NP ) 

CALL  HP (SF, SPF, SWEEP) 

IDELY  =  14 

DO  42  ID  =  1, IDELY 
CALL  KEITH (MV) 
42    CONTINUE 

DO  40  J=1,NP 
AVG=0 . 0 
DO  41  JA=1,NA 
300  CALL  KEITH(MV) 

IF(MV. GT.20000.OR.MV.LT. 100)  GOTO  300 
AVG=AVG+FLOAT ( MV ) 
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4 1         CONTINUE 

AMP=AVG / FLOAT ( NA ) 
DY=AMP 

DX=SLOPE*FLOAT(J)+SF 
WRITE (50,*)  DX,DY 
C  IF(AMP.GT.YMAX)  THEN 

C  YMAX=AMP 

C  IMAX=J 

C  ENDIF 

C  IF  (AMP.LT.YMIN)  YMIN=AMP 

40    CONTINUE 
C       CALL  PLOT(NP,DX,DY,XMAX,XMIN,YMAX,YMIN) 
CALL  KEITHA(T) 
CALL  PRES(P) 
WRITE (50,*)  T,P 
IG=IG-1 
CLOSE (50) 
return 
End 


PROGRAM  GETSPG 

$ LARGE 

IMPLICIT  REAL*8  (A-H,0-Z) 

REAL*8  PA(10000) ,TA(10000) ,FR(10000) ,CA,DSQRT 

REAL*8  B(100) ,B1(100) ,B2(100) ,T(100) 

REAL*8  C(100) ,C1(100) ,C2 ( 100 ) ,BD( 100 ) ,N 

INTEGER  TIM 

ZERO=0.0D  00 

ONE=1.0D  00 

TWO=2.0D  00 

EPK=119.8D  00 

B0=49.80D-03 

CON=6.8046D-02 

KEL=273.15 

OPEN  ( 10 , FILE= 'ARGON . DAT ' , STATUS= ' OLD ■ ) 

1=1 
200      READ(10,*,END=100)  Tl ,T2 ,Pl ,P2 , AMP,FR( I ) ,TIM, IN 

IF(Pl.EQ.O.ODO)  GOTO  200 

IF(P2.EQ.0.0D0)  GOTO  200 

IF(Tl.EQ.O.ODO)  GOTO  200 

IF(T2.EQ.0.0D0)  GOTO  200 

TA(I)=(Tl+T2)/TWO 

PA ( I ) = ( P1+P2 ) *CON/TWO 

1=1+1 

GOTO  200 
100      CLOSE (10) 

M=I-1 

OPEN  ( 10 , FILE= ' TVC . DAT ' , STATUS= ' OLD ' ) 

DO  25  1=1,74 

READ(10,*)  T(I),B(I),B1(I),B2(I), 
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CBD(I),C(I),C1(I),C2(I) 
25      CONTINUE 

CLOSE  (10) 

OPEN  (30,FILE='NAT.GAS' , STATUS= ' OLD • ) 
OPEN  (40, FILE= ' SPEED . GAS ' , STATUS= ' NEW ' ) 
c         open  ( 50, f ile='a.gas ' , status= ' new' ) 
c         open  (60, f ile='b.gas ' , status='new' ) 
c         open  (70, f ile=*c.gas ' , status=' new' ) 
N=ZERO 
500     READ(30,*,END=300)  T1,T2 ,P1,P2,AMP,FREQ,TIM, IN 
c         write  (60,*)  amp, in 
c         write  (50,*)  freq 
N=N+ONE 

IF(Pl.EQ.0.0.OR.P2.EQ.0.0)  GOTO  500 
TEMP=(T1+T2)/TW0 
P=(P1+P2)/TW0 
DT=T2-T1 
DP=P2-P1 

IF(TA(1) .LT.TEMP)  GOTO  500 
IF(TA(M) .GT.TEMP)  GOTO  500 
IF(AMP.GT. 19000. 0D0)  IN=IN-1 
IF(AMP.LT.7000.0D0)  IN=IN+1 
c         write  (60,*)  freq 
CALL  AMPGET(AMP,IN) 
c         write  (70,*)  freq 
DO  20  1=1, M 
ISTOP=0 

IF  (TA(I) .LT.TEMP)  GOTO  400 
ISTOP=l 
20      CONTINUE 
400      IF  (ISTOP.EQ.l)  GOTO  500 

FAR=(FR(I)-FR(I-1) ) / (TA( I ) -TA( 1-1 ) )* 
C(TEMP-TA(I-1) )+FR(I-l) 

PAR=(PA(I)-PA(I-1) )/(TA(I)-TA(I-l) )* 
C(TEMP-TA(I-1) )+PA(I-l) 
TEMP=TEMP+KEL 
TS=TEMP/EPK 

CALL  QAND(BS,BS1,BS2,BSD,CS,CS1, 
CCS2,TS,T,B,B1,B2,BD,C,C1,C2) 
C         WRITE (*,*)  '  QAND' 
BV=BS*B0 
CV=CS*B0*B0 

CALL  VERVOL( PAR, V, TEMP, BV,CV) 
C  WRITE (*,*)  '  VERVOL* 

VS=V/B0 

CALL  SPEED ( CAR , TEMP , BS , BS 1 , BS2 , BSD , CS , CS 1 , CS2 , VS ) 
C  WRITE (*,*)  '  SPEED' 

CG=FREQ*CAR/FAR 
TEMP=TEMP-KEL 

WRITE (40,*)  N, DP, DT,P, TEMP, AMP, FREQ, CG, TIM, CAR, FAR 
GOTO  500 
close(70) 
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close(60) 
close(50) 
300  CLOSE (40) 
CLOSE (30) 
END 

SUBROUTINE  QAND(BS, BS1 , BS2 ,BSD, 
CCS,CS1,CS2,TS,T,B,B1,B2,BD,C,C1,C2) 
IMPLICIT  REAL*8  (A-H,0-Z) 
REAL*8  6(1) ,61(1),62(1),T(1) ,0(1) ,C1(1) ,02(1) ,6D(1) 
REAL*8  M 
DO  20  1=2,74 
IF  (TS.GT.T(I-l) .AND.TS.LT.T(I) )  THEN 
M=(TS-T(I-1) )/(T(I)-T(I-l) ) 
6S=6(I-1)+M* (6(1) -6(1-1) ) 
8S1=81(I-1)+M*(81(I)-81(I-1) ) 
8S2=82(I-1)+M*(82(I)-82(I-1) ) 
BSD=8D(I-1)+M*(BD(I)-8D(I-1) ) 
CS=C(I-1)+M*(C(I)-C(I-1) ) 
CS1=C1(I-1)+M*(C1(I)-C1(I-1) ) 
CS2=C2(I-1)+M*(C2(I)-C2(I-1) ) 
RETURN 
ENDIF 
20      CONTINUE 

WRITE (*,*) 'TSTAR  OUT  OF  RANGE' 
RETURN 
END 
C 

SUBROUTINE  VERVOL(P, V,T,B,CV) 
IMPLICIT  REAL*8  (A-H,0-Z) 
R=8.20575D-02 
TOL=1.0D-16 
C  INPUT  P  IN  ATM 
C  INPUT  T  IN  KELVIN 
V=R*T/P 
10     VN=R*T/P*(1.0D  00  +  B/V  +  CV/V/V) 
TEST=V/VN 

IF ( TEST. GT. 1.0)  THEN 
TEST=1.0D  00  -  1.0D  00/TEST 
ELSE 

TEST=1.0D  00  -TEST 
ENDIF 
V=VN 

IF(TEST.GT.TOL)   GOTO  10 
RETURN 
END 
C 

SUBROUTINE  SPEED ( C , T , BS , BS 1 , BS2 , BSD , OS , CS 1 , CS2 , VS ) 

IMPLICIT  REAL*8  (A-H,0-Z) 

REAL* 8  M,DSQRT 

M=39.948D-03 

R=8.31441D  00 
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GAMA=5 . 0D0/2 . 0D0-BS2 /VS+ 
C(BSD*BSD-CS+CS1-0.5D0*CS2)/VS/VS 

GAMA=GAMA/ ( 3 . 0D0/2 . ODO- ( 2 . 0D0*BS1+BS2 ) /VS- 
C ( 2 . 0D0*CS1+CS2 ) /2 . ODO/VS/VS ) 

C=GAMA*R*T/M* ( 1 . OD  00+2 . 0D0*BS/VS+3 . ODO*CS/VS/VS ) 

C=DSQRT(C) 

RETURN 

END 


PROGRAM  TOTDAT 


ARGC,ARGF 
X6) 


CLOSE 
CLOSE 
END 


ARGF , IT , ARGC , EXPC , THEOC 
,2F9.2,I8,3F9.3) 


COMMON  GROUP  1 

integer 
integer 
integer 
integer 
integer 
integer 
integer 
integer 
integer 
integer 
integer 
integer 

COMMON  GROUP  2 


UNL,  UNT,  GTL,  SDC,  PPC,  GET,  TCT 
LLO,  DCL,  PPU,  SPE,  SPD,  PPE,  PPD 
ERR,  TIMO,  END,  SRQI,  RQS,   CMPL,  LOK 
REM,  CIC,   ATN,  TACS,  LACS,  DTAS ,  DCAS 
EDVR,  ECIC,  ENOL,  EADR,  EARG,  ESAC,  EABO 
ENEB,  EOIP,  ECAP,  EFSO,  EBUS,  ESTB,  ESRQ 
BIN,  XEOS,  REOS 
TNONE,   TlOus,  T30us,  TlOOus,  T300us 
Tims,    T3ms,   TlOms,  T30ms,   TlOOms 
T300ms,  Tls,    T3s,    TlOs,    T30s 
TlOOs,   T300s,  TlOOOs 
S,LF 


character*50  flname 

integer*4      bd,dvm,v,cnt ,mask 

integer*4      spr,ppr 

data  UNL/63/,UNT/95/,GTL/01/,SDC/04/,PPC/05/ 

data  GET/08/, TCT/09/, LLO/ 17 /,DCL/20/,PPU/21/ 

data  SPE/24/,SPD/25/,PPE/96/,PPD/112/ 

data  TIMO/ 16384/, END/8 192 /,SRQI/4096/ 

data  RQS/2048/,CMPL/256/,LOK/128/, REM/64 /,CIC/32/ 
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data  ATN/ 16/, TACS/8/, LACS/4 / ,DTAS/2/,DCAS/l/ 

data  EDVR/0/,ECIC/l/,ENOL/2/,EADR/3/,EARG/4/ 

data  ESAC/5/,EABO/6/,ENEB/7/,EOIP/10/,ECAP/ll/ 

data  EFSO/ 1 2 / , EBUS / 14 / , ESTB/ 1 5 / , ESRQ/ 16/ 

data  BIN/4096 /,XE0S/2 048/ ,REOS/ 1024/ 

data  TNONE/0/,T10us/1/,T30us/2/,T100us/3/,T300us/4/ 

data  Tlms/5/,T3ms/6/fTl0ms/7/,T30ms/8/,Tl00ms/9/ 

data  T300ms/10/,Tls/ll/,T3s/12/Tl0s/13/,T30s/14/ 

data  T100s/15/,T300s/16/,T1000s/17/ 

data  S/08/,LF/10/ 
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